Datahack Logo Logo de Python

Profesor: Julio Antonio Soto

Scikit-Learn: machine learning hecho sencillo con Python

Ahora que ya tenemos una base del lenguaje Python y conocemos las posibilidades de Numpy y pandas (así como de visualización), ha llegado la hora de aplicar todo lo aprendido al mundo del aprendizaje automático.

Scikit-Learn es, probablemente, la biblioteca de machine learning más extensa que hay. Rápida, estable y en continua evolución, Scikit-Learn (se abrevia como sklearn) es la respuesta evidente a hacer machine learning con Python.

Sklearn está escrita en una mezcla de Python, Cython, C y C++. Utiliza extensamente Numpy y Scipy por debajo, y su integración con nuestras Series y DataFrames de pandas es tan sencilla como te imaginas.

Índice

  1. Introducción
  2. Cosas importantes antes de ponerse a machacar datos
    1. Creación de datos sintéticos
    2. Datasets precargados
    3. Separación de datos de training y testing
    4. Feature standarization
  3. Clasificación
    1. Preparemos un dataset sencillo e intuitivo
    2. Funcionamiento de los clasificadores
    3. Separamos el dataset en training y testing
    4. Entrenamiento de clasificadores
      1. Regresión logística
      2. Árboles
      3. Nearest Neighbors
      4. Naïve Bayes
      5. Otros clasificadores: SVMs, redes y demás
  4. Ensemble Learning
    1. Bagging
    2. Random Forests
    3. AdaBoost
    4. Gradient Boosting
  5. Métricas de clasificación
    1. Confusion matrix
    2. Precision, recall y F1-score
      1. Precision
      2. Recall
      3. F1-score
    3. Curva ROC
  6. Regresión
    1. Regresión lineal simple (mínimos cuadrados ordinarios)
    2. Regularización: Ridge, LASSO y Elastic Net
    3. Regresión polinómica: transformación no lineal de nuestro dataset
    4. Otros regresores
      1. Árboles regresores
      2. Ensembles
      3. Regresión por Nearest Neighbors
      4. Support Vector Regression
  7. Clustering
    1. Los datos
    2. Algoritmos de clustering
      1. K-Means
      2. Gaussian Mixture Model
      3. Affinity Propagation
      4. Mean Shift
      5. DBSCAN
  8. Feature Engineering
    1. Reducción de dimensionalidad
      1. Principal Component Analysis
    2. Selección de variables
      1. Filtro de variables por varianza
      2. Selección de variables univariante
      3. Uso de árboles u otros clasificadores y regresores para seleccionar variables
      4. Recursive Feature Elimination
  9. Validación y selección de modelos
    1. Training, validación, testing
    2. Cross Validation
      1. Validación cruzada empotrada en la especificación de los modelos
      2. GridSearchCV
  10. Pipelines
  11. Bibliografía y recursos adicionales

1. Introducción

Sklearn es una biblioteca que contiene una gran variedad de modelos escritos, y listos para ser entrenados y utilizados:

  • Modelos de clasificación
  • Modelos de regresión
  • Modelos de clustering
  • Modelos de reducción de dimensionalidad
  • Utilidades para validación y selección de modelos
  • Utilidades para preprocesado de datos y feature engineering
  • Una arquitectura de pipelines para hacer workflows complejos de todos los puntos anteriores

A todo el mundo que conozco le encanta el siguiente diagramita que aparece en la página web de sklearn, sobre "cómo elegir tu modelo ideal":

mapa de sklearn

Evidentemente se trata de una ultra-generalización, pero no por ello deja de ser bonita. Y no se puede decir que sea errónea tampoco...

Iremos haciendo referencias a la documentación de sklearn. Solo diré que es probablemente la mejor documentación que he visto, y solo con ella puedes aprender muchísimo sobre machine learning. No, no es broma.

Suficiente. Pongámonos manos a la obra.

In [1]:
%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
plt.style.use("ggplot")

2. Cosas importantes antes de ponerse a machacar datos

Estarás deseando empezar a lanzar modelos... Pero antes vamos a ver algunas de las opciones que nos ayudarán a hacerlo:

2.A. Creación de datos sintéticos

A veces... no tenemos datos. Es lo que hay. Pero eso no tiene que pararnos los pies, sobre todo cuando se trata de aprender. De hecho, la mayoría de I+D en machine learning e inteligencia artificial se hace con datos sintéticos (es decir, generados de forma aleatoria por nosotros); o bien con datasets muy conocidos.

Afortunadamente, scikit-learn trae de serie una serie de funciones que nos permiten generar datos sintéticos de forma muy sencilla.

Por ejemplo: la función make_blobs genera muestras gausianas (es decir, distribuciones normales en dos o más dimensiones) isotrópicamente distribuidas. Vamos a crear tres blobs:

In [2]:
# La biblioteca scikit-learn se llama sklearn
# (para el tema de imports y tal):
from sklearn.datasets import make_blobs

mi_primer_dataset = make_blobs(n_samples=10,     # Numero de puntos que creamos
                               n_features=2,     # Numero de atributos (variables)
                               centers=3,        # Numero de centros (numero de gausianos) -> 3 blobs
                               cluster_std=1.5,  # Desviacion tipica de los datos
                               random_state=2)   # Si queremos poner un "seed" al generador de números aleatorios

mi_primer_dataset
Out[2]:
(array([[  0.96455381,   0.46894968],
        [ -1.21779287, -11.15836353],
        [  1.80183704,  -2.1877917 ],
        [ -2.10087691,  -3.74757963],
        [ -0.45292089,  -6.04316334],
        [ -0.52577983, -11.34940749],
        [ -0.12855687,  -1.28001427],
        [ -1.20778828,  -4.87647215],
        [ -2.9098058 ,  -3.62795484],
        [ -2.86703029, -10.84498679]]), array([1, 0, 1, 2, 0, 0, 1, 2, 2, 0]))

Nos devuelve una tupla con dos arrays de Numpy: la primera son las coordenadas de cada punto, y la segunda a qué centro pertenece cada punto; de forma que el punto [ 0.96455381, 0.46894968] pertenece al centro 1, el punto [ -1.21779287, -11.15836353] pertenece al centro 0, y así.

Dibujemos con matplotlib. Preparamos un DataFrame sencillito con los datos (para verlo más claro; aunque podríamos dibujar a pelo con matplotlib y las arrays de Numpy):

In [3]:
# Primero vamos a meter en un DataFrame las coordenadas de cada punto:
dataset_blobs = pd.DataFrame(mi_primer_dataset[0]) # Primer elemento de la tupla

# Y vamos a renombrar las columnas:
dataset_blobs.columns = ["x", "y"]

# Y podemos añadir el centro al que pertenece cada uno al DataFrame.
# ¿Por qué no? Es común llamarlo "label":
dataset_blobs["label"] = mi_primer_dataset[1]

dataset_blobs
Out[3]:
x y label
0 0.964554 0.468950 1
1 -1.217793 -11.158364 0
2 1.801837 -2.187792 1
3 -2.100877 -3.747580 2
4 -0.452921 -6.043163 0
5 -0.525780 -11.349407 0
6 -0.128557 -1.280014 1
7 -1.207788 -4.876472 2
8 -2.909806 -3.627955 2
9 -2.867030 -10.844987 0
In [4]:
# Los colores por defecto de Matplotlib cuando
# hay distintos puntos son feos...
import matplotlib.cm as color_maps
# Ahora podemos usar estos:
# http://matplotlib.org/examples/color/colormaps_reference.html#color-example-code-colormaps-reference-py

dataset_blobs.plot(kind="scatter",
                   x="x",
                   y="y",
                   c="label",
                   cmap=color_maps.brg,
                   figsize=(8,6))
pass

Sosísimo. Muy pocos puntos. Vamos a hacer lo mismo, pero aumentando el número de puntos:

In [5]:
mi_primer_dataset = make_blobs(n_samples=250, # 250 puntos será más vistoso.  
                               n_features=2,    
                               centers=3,       
                               cluster_std=0.5,
                               random_state=2)

dataset_blobs = pd.DataFrame(mi_primer_dataset[0])
dataset_blobs.columns = ["x", "y"]
dataset_blobs["label"] = mi_primer_dataset[1]
dataset_blobs.plot(kind="scatter",
                   x="x",
                   y="y",
                   c="label",
                   cmap=color_maps.brg,
                   figsize=(8,6))
pass

Mejor. Hacer un clustering con este dataset sería coser y cantar.

El resto de generadores de datos sintéticos funcionan prácticamente igual que make_blobs. Aquí está la lista completa.

2.B. Datasets precargados

Pero siempre podemos preferir tirar de un dataset real. Si no tenemos el nuestro propio, sklearn viene con unos cuantos precargados, así como con utilidades para descargar automágicamente de repositorios online como puede ser mldata.org.

Vamos a cargar un famoso dataset de precios de casas en Boston:

In [6]:
from sklearn.datasets import load_boston

datos_boston = load_boston()

Con load_boston() hemos cargado un diccionario de Python, que tiene las siguientes claves:

In [7]:
datos_boston.keys()
Out[7]:
dict_keys(['DESCR', 'feature_names', 'target', 'data'])

En DESCR supongo que estará la descripción del dataset...

In [8]:
print(datos_boston["DESCR"])
Boston House Prices dataset
===========================

Notes
------
Data Set Characteristics:  

    :Number of Instances: 506 

    :Number of Attributes: 13 numeric/categorical predictive
    
    :Median Value (attribute 14) is usually the target

    :Attribute Information (in order):
        - CRIM     per capita crime rate by town
        - ZN       proportion of residential land zoned for lots over 25,000 sq.ft.
        - INDUS    proportion of non-retail business acres per town
        - CHAS     Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
        - NOX      nitric oxides concentration (parts per 10 million)
        - RM       average number of rooms per dwelling
        - AGE      proportion of owner-occupied units built prior to 1940
        - DIS      weighted distances to five Boston employment centres
        - RAD      index of accessibility to radial highways
        - TAX      full-value property-tax rate per $10,000
        - PTRATIO  pupil-teacher ratio by town
        - B        1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
        - LSTAT    % lower status of the population
        - MEDV     Median value of owner-occupied homes in $1000's

    :Missing Attribute Values: None

    :Creator: Harrison, D. and Rubinfeld, D.L.

This is a copy of UCI ML housing dataset.
http://archive.ics.uci.edu/ml/datasets/Housing


This dataset was taken from the StatLib library which is maintained at Carnegie Mellon University.

The Boston house-price data of Harrison, D. and Rubinfeld, D.L. 'Hedonic
prices and the demand for clean air', J. Environ. Economics & Management,
vol.5, 81-102, 1978.   Used in Belsley, Kuh & Welsch, 'Regression diagnostics
...', Wiley, 1980.   N.B. Various transformations are used in the table on
pages 244-261 of the latter.

The Boston house-price data has been used in many machine learning papers that address regression
problems.   
     
**References**

   - Belsley, Kuh & Welsch, 'Regression diagnostics: Identifying Influential Data and Sources of Collinearity', Wiley, 1980. 244-261.
   - Quinlan,R. (1993). Combining Instance-Based and Model-Based Learning. In Proceedings on the Tenth International Conference of Machine Learning, 236-243, University of Massachusetts, Amherst. Morgan Kaufmann.
   - many more! (see http://archive.ics.uci.edu/ml/datasets/Housing)

Pues sí. feature_names nos da los nombres de las columnas, data son las variables independientes, y target es el valor dependiente (el que queremos predecir, es decir: el valor de cada casa). Cargarlo todo en un DataFrame es coser y cantar:

In [9]:
# Cargamos las columnas en un DataFrame:
dataset_boston = pd.DataFrame(datos_boston["data"])

# Y ponemos los nombres a las columnas
# (damos por hecho que el orden de las columnas
# es el correcto):
dataset_boston.columns = datos_boston["feature_names"]

# Y añadimos la columna de precios
# (que si leemos más arriba en la descripción,
# podemos ver que quieren que la columna se llame
# MEDV. Pero somos unos rebeldes):
dataset_boston["Precio"] = datos_boston["target"]

dataset_boston[:10]
Out[9]:
CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX PTRATIO B LSTAT Precio
0 0.00632 18.0 2.31 0.0 0.538 6.575 65.2 4.0900 1.0 296.0 15.3 396.90 4.98 24.0
1 0.02731 0.0 7.07 0.0 0.469 6.421 78.9 4.9671 2.0 242.0 17.8 396.90 9.14 21.6
2 0.02729 0.0 7.07 0.0 0.469 7.185 61.1 4.9671 2.0 242.0 17.8 392.83 4.03 34.7
3 0.03237 0.0 2.18 0.0 0.458 6.998 45.8 6.0622 3.0 222.0 18.7 394.63 2.94 33.4
4 0.06905 0.0 2.18 0.0 0.458 7.147 54.2 6.0622 3.0 222.0 18.7 396.90 5.33 36.2
5 0.02985 0.0 2.18 0.0 0.458 6.430 58.7 6.0622 3.0 222.0 18.7 394.12 5.21 28.7
6 0.08829 12.5 7.87 0.0 0.524 6.012 66.6 5.5605 5.0 311.0 15.2 395.60 12.43 22.9
7 0.14455 12.5 7.87 0.0 0.524 6.172 96.1 5.9505 5.0 311.0 15.2 396.90 19.15 27.1
8 0.21124 12.5 7.87 0.0 0.524 5.631 100.0 6.0821 5.0 311.0 15.2 386.63 29.93 16.5
9 0.17004 12.5 7.87 0.0 0.524 6.004 85.9 6.5921 5.0 311.0 15.2 386.71 17.10 18.9

Distintos datasets precargados pueden venir en distintos formatos; pero siempre utilizando una estructura de dato de Python (probablemente algún tipo de colección) que ya conocemos, así que no deberíamos tener problemas en cargar cualquiera de ellos.

2.C. Separación de datos de training y testing

Es común no utilizar todos los datos de los que disponemos para entrenar nuestros modelos. La causa es bien sencilla: los algoritmos de aprendizaje automático tienden a "aprenderse de memoria los datos" con los que entrenan. Cuando entrenamos un modelo, lo más frecuente es que dicho modelo intente minimizar el error cometido dentro de ese dataset de entrenamiento ($E_{in}$, de "dentro" de los datos mostrados al modelo). La gracia del aprendizaje automático es que nuestros modelos sean capaces de generalizar, es decir, que sean capaces de predecir bien en datos a los que no han sido expuestos anteriormente (el error cometido con datos "nuevos" con los que un modelo no ha entrenado se puede escribir como $E_{out}$).

En un mundo ideal los modelos generalizan perfectamente, de forma que $E_{in} \approx E_{out}$, es decir: el error cometido con los datos de entrenamiento es (a ser posible) pequeño, y es igual de pequeño con datos nuevos.

Verificar la diferencia entre $E_{in}$ y $E_{out}$ es muy importante para ver la precisión final que tendrá nuestro modelo. Por eso es buena costumbre reservar parte de los datos y apartarlos del entrenamiento. De forma que:

  • Utilizaremos una parte de los datos llamada training set para entrenar los modelos.
  • Utilizaremos el resto de los datos (test set) para ver cómo de bueno es nuestro modelo con datos con los que no ha entrenado.

Más adelante entraremos en el tema de la validación.

El caso es que sklearn nos permite dividir nuestro dataset en training y test sets de una forma muy sencilla, con sklearn.model_selection.train_test_split.

Ya que tenemos el dataset de los precios de Boston cargado, vamos a separarlo en training set (80% de los datos) y test set (20% de los datos):

In [10]:
from sklearn.model_selection import train_test_split

datos_spliteados = train_test_split(dataset_boston,
                                    train_size=0.8, # 80% training
                                    test_size=0.2   # 20% testing
                                   )

datos_spliteados es ahora una tupla que contiene los dos DataFrames: el de training y el de testing:

In [11]:
boston_training_set = datos_spliteados[0]
boston_test_set = datos_spliteados[1]

print("Training set (filas, columnas): " + str(boston_training_set.shape))
print("Test set (filas, columnas): " + str(boston_test_set.shape))
Training set (filas, columnas): (404, 14)
Test set (filas, columnas): (102, 14)

2.D. Feature standarization

Es frecuente encontrarse con datasets que tienen variables que pertenecen a dominios completamente distintos, y con variopintas magnitudes. Por muy inteligentes que sean nuestros modelos, mezclar peras con manzanas no siempre es ideal. Por ejemplo: si queremos hacer un modelo de scoring de riesgo para clientes, nos podemos encontrar con variables como edad (que probablemente tome valores entro 18 y 100 años) y nivel de ingresos (que estará en euros/miles de euros). Distintas magnitudes y, sobre todo, distintos valores absolutos. En un caso como este es posible que los modelos decidan "dar más importancia" a la variable ingresos; que seguro va a tener unos valores máximos y mínimos más extremos, y probablemente más varianza. Seguro que también va a haber outliers más drásticos...

Una buena solución para esto es lo que se llama feature standarization. Feature standarization es una técnica que sirve para que todas las variables hablen el mismo idioma, y tomen valores relativos similares entre ellas. Además, por lo general, ayuda a que los modelos converjan (o terminen siendo precisos) en menos tiempo e iteraciones.

Hay varias formas de hacer feature standarization. No obstante, la más utilizada es la que hace que todas las variables tengan media $0$ y desviación típica $1$. Y resulta que sklearn tiene una función que hace exactamente esto. En la bibliografía e Internet podrás encontrar sitios donde ponen que: feature standarization = feature demeaning + feature scaling (es decir, hacer dos cosas: quitar media y escalar varianza/desv.típica).

sklearn.preprocessing incluye varias funciones para hacer este tipo de transformaciones iniciales a nuestras variables. Para conseguir que las variables tengan media $0$ y desviación típica $1$ utilizamos sklearn.preprocessing.StandardScaler.

Continuemos con los datos de Boston: vamos a coger el training set, y ver la media y desviación típica de cada variable. Con un .describe() en el DataFrame bastará:

In [12]:
boston_test_set.describe()
Out[12]:
CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX PTRATIO B LSTAT Precio
count 102.000000 102.000000 102.000000 102.000000 102.000000 102.000000 102.000000 102.000000 102.000000 102.000000 102.000000 102.000000 102.000000 102.000000
mean 3.482922 13.696078 11.443725 0.088235 0.553215 6.372412 67.395098 3.706814 9.568627 411.745098 18.008824 347.532157 12.516667 23.266667
std 8.046422 26.461737 7.043995 0.285037 0.121537 0.687400 29.700093 2.048705 8.833453 171.868686 2.268646 104.684729 7.679847 10.768282
min 0.009060 0.000000 0.460000 0.000000 0.394000 4.906000 6.600000 1.174200 1.000000 188.000000 13.000000 24.650000 1.730000 5.000000
25% 0.085432 0.000000 5.177500 0.000000 0.445000 5.935000 41.950000 1.990325 4.000000 276.000000 16.100000 372.247500 6.292500 15.725000
50% 0.257780 0.000000 10.010000 0.000000 0.545500 6.261000 76.450000 2.970350 5.000000 332.500000 18.350000 390.590000 11.360000 21.700000
75% 3.245900 20.000000 18.100000 0.000000 0.629250 6.798500 94.725000 5.438450 24.000000 666.000000 20.200000 396.750000 17.192500 28.675000
max 67.920800 100.000000 27.740000 1.000000 0.871000 8.297000 100.000000 8.792100 24.000000 711.000000 21.200000 396.900000 34.770000 50.000000

Podemos ver que tenemos variopintas medias y desviaciones típicas en cada variable. Vamos a estandarizarlas con StandardScaler:

In [13]:
from sklearn.preprocessing import StandardScaler

# Creamos una instancia de StandardScaler:
scaler = StandardScaler()

# Aplicamos el método .fit(), pasando como
# argumento nuestro dataset:
scaler.fit(boston_training_set)

# Ese .fit() computa la media y desviación
# típica de cada variable, dejando preparado
# el escalador. Para aplicarlo sobre nuestros
# datos, llamamos ahora al método .transform()
# (ha de ser después de .fit(), claro):
boston_training_set_scaled = scaler.transform(boston_training_set)

boston_training_set_scaled debería tener ahora nuestras mismas variables, pero con media $0$ y desviación típica $1$. Vamos a verlo:

In [14]:
# Resulta que scaler.transform() devuelve un array/matriz de Numpy,
# no un DataFrame de pandas. Podemos coger las columnas del
# boston_training_set original, y poner los nombres:
boston_training_set_scaled = pd.DataFrame(boston_training_set_scaled)
boston_training_set_scaled.columns = boston_training_set.columns

boston_training_set_scaled.describe()
Out[14]:
CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX PTRATIO B LSTAT Precio
count 4.040000e+02 4.040000e+02 4.040000e+02 4.040000e+02 4.040000e+02 4.040000e+02 4.040000e+02 4.040000e+02 4.040000e+02 4.040000e+02 4.040000e+02 4.040000e+02 4.040000e+02 4.040000e+02
mean 4.080894e-17 -7.763317e-17 2.836015e-16 -1.096483e-16 4.704707e-16 1.201459e-15 -4.554937e-17 -5.166384e-17 -6.128211e-17 1.385031e-16 9.057661e-16 -9.868344e-16 1.511442e-18 -5.867144e-17
std 1.001240e+00 1.001240e+00 1.001240e+00 1.001240e+00 1.001240e+00 1.001240e+00 1.001240e+00 1.001240e+00 1.001240e+00 1.001240e+00 1.001240e+00 1.001240e+00 1.001240e+00 1.001240e+00
min -4.141983e-01 -4.803533e-01 -1.514983e+00 -2.622653e-01 -1.486405e+00 -3.833815e+00 -2.378291e+00 -1.268306e+00 -9.849157e-01 -1.314101e+00 -2.810699e+00 -4.100056e+00 -1.459754e+00 -1.982424e+00
25% -4.057514e-01 -4.803533e-01 -8.616747e-01 -2.622653e-01 -8.483833e-01 -5.385316e-01 -8.317638e-01 -7.958263e-01 -6.391110e-01 -7.535172e-01 -5.502032e-01 1.960550e-01 -7.843453e-01 -5.911008e-01
50% -3.855532e-01 -4.803533e-01 -3.669219e-01 -2.622653e-01 -1.491818e-01 -9.078743e-02 3.110079e-01 -2.727843e-01 -5.238427e-01 -4.612979e-01 2.503891e-01 3.728721e-01 -1.889429e-01 -1.597050e-01
75% 6.079083e-03 7.691391e-02 1.033654e+00 -2.622653e-01 6.024599e-01 4.697797e-01 9.022211e-01 6.131643e-01 1.666254e+00 1.542492e+00 7.684193e-01 4.249713e-01 6.032601e-01 3.031170e-01
max 9.778565e+00 3.754878e+00 2.448911e+00 3.812933e+00 2.761245e+00 3.572767e+00 1.122123e+00 3.921014e+00 1.666254e+00 1.810856e+00 1.616105e+00 4.334592e-01 3.612060e+00 3.160043e+00

Las medias son ahora valores muy cercanos a $0$; y la desviaciones típicas valores muy cercanos a $1$. Por tanto, nuestras variables ahora están estandarizadas.

Nota: con el fin de simplificar el código, en este ejemplo hemos estandarizado también la target variable, variable dependiente o variable a predecir (o sea, el precio de las casas). En un escenario real no vas a querer hacer esto; solo vas a querer estandarizar las variables que no quieres predecir.

Sklearn incluye otros varios transformadores y preprocesadores de datos. Puedes leer más sobre ellos aquí y aquí.

Con conocimientos de Numpy, pandas y sklearn.preprocessing no deberíamos tener problemas para preparar nuestros datasets para el machine learning que va luego.

Sin más dilación, empecemos con los modelos...

3. Clasificación

Desde mi experiencia, diría que el 80% de los problemas de machine learning/data science que nos encontramos por ahí son problemas de clasificación. Los seres humanos tendemos a "etiquetar" todo dentro de una serie de grupos controlados, con el fin de hacer las cosas más sencillas para nuestro consciente. De hecho, nos pasamos el día clasificando cosas de forma involuntaria. Gracias a eso, podemos ver dos mesas distintas; pero sabemos que ambas son mesas (clasificamos las dos como mesas; diferenciándolas de otras labels: sillas, camas...).

Los problemas de clasificación son sencillos de entender: tenemos un dataset con una serie de variables independientes o features, y una variable que es la clasificación, clase o label de cada observación. Con esos datos, entrenamos en modelo. Cuando en el futuro nos lleguen únicamente las features, seremos capaces de predecir la label de cada observación.

Los problemas de clasificación se resuelven mediante apendizaje supervisado. Esto quiere decir que tenemos en nuestro dataset de entrenamiento las labels de cada observación, y las utilizamos para "enseñar" al modelo.

3.A. Preparemos un dataset sencillo e intuitivo

Vamos a crear un dataset simple para recrear un problema de clasificación en dos dimensiones (¡que podamos dibujar!). Vamos a generar datos sintéticos. Podríamos utilizar make_blobs para este propósito. Pero los desarrolladores de sklearn son muy listos, y saben que la clasificación es el problema más común en machine learning. De forma que han creado una función que genera automáticamente datasets sintéticos para problemas de clasificación: sklearn.datasets.make_classification

In [15]:
from sklearn.datasets import make_classification

dataset = make_classification(n_samples=600, # Número de observaciones
                              n_features=2, # Número de features
                              n_redundant=0, # Número de features "redundantes" que añaden ruido
                              n_classes=2, # Número de distintas clases o distintas labels
                              random_state=749 # Seed del generador de números aleatorios
                             )

Al igual que con make_blobs, el resultado es una tupla con dos arrays de Numpy: la primera tiene las features, y la segunda las labels:

In [16]:
# Metemos en un DataFrame las features primero:
dataset_clasificacion = pd.DataFrame(dataset[0])

# Y ahora las labels:
dataset_clasificacion["label"] = dataset[1]

dataset_clasificacion[:5]
Out[16]:
0 1 label
0 0.574208 1.609926 1
1 -0.666932 -0.698694 0
2 0.452082 -2.546917 1
3 -0.854633 -0.364912 0
4 -0.277625 -0.411260 0

Ahora inventémonos una historia para nuestros datos. Supongamos que este dataset muestra en cada fila la observación de los niveles de concentración de dos hormonas masculinas (hormona_a y hormona_b) de un hombre. Y supongamos que la label que acompaña a las medidas de hormonas de cada hombre nos dice si dicho hombre se ha quedado calvo o no (donde $0$ es hombre con pelo, y $1$ es hombre calvo). Vamos a intentar, por tanto, predecir si un hombre se quedará calvo o no a partir de sus niveles hormonales.

Disclaimer: DataHack no se hace responsable del posible mal gusto del ejemplo de clasificación elegido.

In [17]:
# Vamos a renombrar las columnas de features:
dataset_clasificacion.columns = ["hormona_a", "hormona_b", "label"]

# Y sumar 4 a los valores de ambas features
# (porque hay concentraciones negativas
# de hormonas si no, y es raro):
dataset_clasificacion["hormona_a"] = dataset_clasificacion["hormona_a"] + 4
dataset_clasificacion["hormona_b"] = dataset_clasificacion["hormona_b"] + 4

dataset_clasificacion[:5]
Out[17]:
hormona_a hormona_b label
0 4.574208 5.609926 1
1 3.333068 3.301306 0
2 4.452082 1.453083 1
3 3.145367 3.635088 0
4 3.722375 3.588740 0
In [18]:
dataset_clasificacion.describe()
Out[18]:
hormona_a hormona_b label
count 600.000000 600.000000 600.000000
mean 4.010365 3.949679 0.500000
std 1.220447 1.267697 0.500417
min 1.033943 1.056018 0.000000
25% 3.074655 3.035652 0.000000
50% 3.708435 3.984949 0.500000
75% 5.002841 4.944744 1.000000
max 8.051309 7.548703 1.000000

Dibujemos:

In [19]:
# Queda mejor si pintamos ceros y unos por separado:
dataset_clasificacion_ceros = dataset_clasificacion[dataset_clasificacion["label"] == 0]
dataset_clasificacion_unos = dataset_clasificacion[dataset_clasificacion["label"] == 1]

axes = dataset_clasificacion_ceros.plot(kind="scatter",
                                        x="hormona_a",
                                        y="hormona_b",
                                        label="No calvo",
                                        c="blue", # Color de los puntos
                                        edgecolors="black", # Color del contorno de los puntos
                                        linewidths=1, # Ancho del contorno de los puntos
                                        s=30, # Tamaño de los puntos
                                        figsize=(8,6))

dataset_clasificacion_unos.plot(kind="scatter",
                                x="hormona_a",
                                y="hormona_b",
                                label="Calvo",
                                c="red",
                                edgecolors="black",
                                linewidths=1,
                                s=30,
                                ax=axes)
pass

Donde cada punto es un señor con sus dos niveles de concentración de hormonas, y si es calvo o no.

3.B. Funcionamiento de los clasificadores

Casi siempre vamos a querer hacer lo mismo:

  1. Cargar un modelo (hacer el import correspondiente)
  2. Entrenarlo con nuestros datos
  3. Probar su rendimiento
  4. Predecir nuevos casos (en este caso, si señores nuevos se quedarán o no calvos)

Pues bien, es muy sencillo. Dado un modelo:

  1. Hacemos from sklearn.loquesea import modelo
  2. modelo.fit(features, labels)
  3. modelo.score(features, labels)
  4. modelo.predict(features, labels)

El orden en el que queremos hacer los pasos puede cambiar, pero por lo general será algo así.

3.C. Separamos el dataset en training y testing

Pues eso. Igual que en el ejemplo de Boston:

In [20]:
dataset_clasificacion_spliteado = train_test_split(dataset_clasificacion,
                                                   train_size=0.8, # 80% training
                                                   test_size=0.2,  # 20% testing
                                                   random_state=22 # seed para el generador de num. al.
                                                  )

dataset_clasificacion_training = dataset_clasificacion_spliteado[0]
dataset_clasificacion_testing = dataset_clasificacion_spliteado[1]

Y dibujamos:

In [21]:
plt.figure(figsize=(16,7))

# Training
dataset_clasificacion_training_ceros = dataset_clasificacion_training[dataset_clasificacion_training["label"] == 0]
dataset_clasificacion_training_unos = dataset_clasificacion_training[dataset_clasificacion_training["label"] == 1]

axes1=plt.subplot(1,2,1, title="Training")

dataset_clasificacion_training_ceros.plot(kind="scatter",
                                        x="hormona_a",
                                        y="hormona_b",
                                        label="No calvo",
                                        c="blue",
                                        edgecolors="black",
                                        linewidths=0.5,
                                        s=30,
                                        ax=axes1)

dataset_clasificacion_training_unos.plot(kind="scatter",
                                x="hormona_a",
                                y="hormona_b",
                                label="Calvo",
                                c="red",
                                edgecolors="black",
                                linewidths=0.5,
                                s=30,
                                ax=axes1)

# Testing
dataset_clasificacion_testing_ceros = dataset_clasificacion_testing[dataset_clasificacion_testing["label"] == 0]
dataset_clasificacion_testing_unos = dataset_clasificacion_testing[dataset_clasificacion_testing["label"] == 1]

axes2 = plt.subplot(1,2,2, title="Testing")
dataset_clasificacion_testing_ceros.plot(kind="scatter",
                                        x="hormona_a",
                                        y="hormona_b",
                                        label="No calvo",
                                        c="blue",
                                        edgecolors="black",
                                        linewidths=0.5,
                                        s=30,
                                        ax=axes2)

dataset_clasificacion_testing_unos.plot(kind="scatter",
                                x="hormona_a",
                                y="hormona_b",
                                label="Calvo",
                                c="red",
                                edgecolors="black",
                                linewidths=0.5,
                                s=30,
                                ax=axes2)
pass

3.D. Entrenamiento de clasificadores

No podemos saber qué modelo de clasificación va a ser el que mejor funcione para un determinado dataset. Podemos tener nuestras hipótesis o corazonadas; pero fiarse solo en eso contradice el método científico...

... De forma que tendremos que probar distintos clasificadores, y ver cuál es el que mejor funciona. Estos son algunos de los clasificadores que incluye sklearn (hay más):

Algoritmo Módulo Multiclase Comentario
Regresión Logística linear_model.LogisticRegresion OvR Emplea la biblioteca liblinear
Árboles de decisión tree.DecisionTreeClassifier Inherente Clasificador de reglas
K Nearest Neighbors neighbors.KNeighborsClassifier Inherente Clasificador en función de la clase de las muestras más cercanas
Naïve Bayes naive_bayes.GaussianNB Inherente Clasificador bayesiano que asume independencia entre atributos
Perceptrón linear_model.perceptron OvR Clasificador linear de margen
SVM (sólo lineal) svm.LinearSVC OvR Emplea la biblioteca liblinear. Sólo kernel lineal
SVM svm.SVC OvO Emplea la biblioteca libsvm. Todo tipo de kernel

Vamos a ir probando entonces:

3.D.a. Regresión logística

A pesar de su nombre la regresión logística no es un algoritmo de regresión, sino de clasificación. Su formulación matemática se basa en la de la regresión lineal, pero luego aplica una función (función logística) para convertir el resultado de dicha regresión en números entre 0 y 1 (por lo general, aunque existen variantes que lo hacen entre -1 y 1). Gracias a este comportamiento podemos interpretar el resultado de predicción como un resultado probabilístico (las probabilidades siempre están entre 0 y 1). Una regla sencilla es que si el resultado > 0.5, se predice como 1, y si el resultado es < 0.5, se predice como 0.

En prácticamente todos los modelos existen una serie de hiperparámetros y argumentos que tenemos que fijar. La formulación matemática de cada modelo será distinta y, por tanto, los resultados también serán distintos, en función de los hiperparámetros que elijamos. La selección de estos hiperparámetros es un mundo complejo, ya que altera la especificación de cada modelo; incluso puede añadir regularizaciones a los mismos. Esta selección de hiperparámetros permite evitar el sobreajuste de los modelos.

Comencemos por entrenar una regresión logística. Evidentemente, con nuestros datos de entrenamiento. El modelo es sklearn.linear_model.LogisticRegression:

In [22]:
from sklearn.linear_model import LogisticRegression

# Creamos una instancia del modelo,
# con todos los hiperparámetros
# y argumentos por defecto:
lr = LogisticRegression()

# Ahora vamos a entrenarla con el conjunto de training:
lr.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]], # nuestras dos features
       y=dataset_clasificacion_training["label"]) # labels
Out[22]:
LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,
          intercept_scaling=1, max_iter=100, multi_class='ovr', n_jobs=1,
          penalty='l2', random_state=None, solver='liblinear', tol=0.0001,
          verbose=0, warm_start=False)

El modelo ha sido entrenado, y vemos que devuelve un montón de opciones e hierparámetros elegidos (que puesto que no hemos elegido nada, ha dejado los que vienen por defecto).

Con el modelo entrenado, vamos a probar a hacer un .predict() sobre los datos de testing, para intentar predecir cuáles están calvos/no calvos. Evidentemente tenemos las labels reales de los datos de testing, y luego ya entraremos en el tema de comprobar rendimiento. De momento queremos ver que con nuevos datos que el modelo no ha visto, nuestra regresión logística predice ceros o unos:

In [23]:
lr.predict(dataset_clasificacion_testing[["hormona_a", "hormona_b"]])
Out[23]:
array([1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0,
       0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0,
       0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1,
       1, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1,
       1, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0,
       1, 0, 0, 0, 0])

Para cada observación de nuestro conjunto de testing, predice 0 o 1. Pero nosotros tenemos los valores reales...

In [24]:
np.array(dataset_clasificacion_testing["label"])
Out[24]:
array([1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0,
       0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 0,
       0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1,
       1, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1,
       1, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0,
       1, 0, 0, 0, 0])

Mirando un poco por encima, vemos que no acierta el 100% en el conjunto de testing. No pasa nada; ahora lo miramos en más profundidad. De hecho, es tan sencillo como utilizar lr.score(features_testing, labels_testing). Básicamente hace lo que hemos hecho nosotros: hace lr.predict() sobre los datos de testing, y compara los resultados de las labels predichas con las reales:

In [25]:
lr.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
         y=dataset_clasificacion_testing["label"])
Out[25]:
0.92500000000000004

Hemos conseguido un número mayor que 0.9. Esto quiere decir que hemos acertado si es 0 o 1 en más del 90% de los casos del dataset de testing (nuestro $E_{out}<10\%$). Esta métrica es conocida como accuracy.

Vale: hemos entrenado la regresión logística y hemos comprobado su accuracy con datos con los que no había entrenado (el conjunto de testing). Si recuerdas la regresión lineal, el resultado del modelo era una función lineal:

$$f(x) = \beta_{0} + \beta_{1}{hormona_a} + \beta_{2}{hormona_b}$$

Pues en la regresión logística es exactamente igual, solo que con la aplicación de la función logística al final. Vamos a ver los parámetros $\beta$ estimados por el modelo:

In [26]:
print(lr.coef_) # beta1 y beta2
print(lr.intercept_) # beta0, también llamado bias o intercept

print("f(x) = " + str(lr.intercept_[0]) + " + " + str(lr.coef_[0][0]) + " hormona_a + " + str(lr.coef_[0][1]) + " hormona_b")
[[ 2.43962443 -0.78176469]]
[-6.29885333]
f(x) = -6.29885332896 + 2.43962442593 hormona_a + -0.781764692684 hormona_b

Esto sería para la regresión lineal normal. Para la regresión logística, tenemos que aplicar la función logística al final (en los libros la veréis también como función sigmoidal o función sigmoide):

$$\theta(f(x)) = \frac{e^{f(x)}}{1+e^{f(x)}}$$

Y así es como conseguimos que nuestro valor predicho $\theta(f(x))$ esté entre 0 y 1.

Te preguntarás si se puede dibujar esto de forma sencilla. La gente de sklearn ha preparado algunas funciones (como esta, esta, esta y esta) para que los clasificadores sean lo más fáciles posible de visualizar. Siéntete libre de utilizar (¡y mejorar!) esta función, que cogí de ellos y alteré ligeramente. Dibuja nuestro dataset, así como las "regiones de decisión" de nuestro clasificador:

In [27]:
def plot_decision_regions(clf, X, y, fig=None, title='', xlabel='', ylabel='', figsize=(8,6)):
    """
    Dibuja las regiones de decisión de un clasificador de sklearn.
    Toma como argumentos:
        - clf: la instancia del clasificador entrenado.
        - X: las features del dataset.
        - y: las labels del dataset.
        - fig: la plt.figure() o axes a usar. Si es None, se creará una figure nueva.
        - title: el título del gráfico.
        - xlabel: el nombre del eje x.
        - ylabel: el nombre del eje y.
        - figsize: si fig es None, el tamaño de la figura a crear (ancho, alto)
    Devuelve:
        Nada. Simplemente dibuja el resultado con Matplotlib.
    """
    from matplotlib.colors import ListedColormap
    cmap_light = ListedColormap(list(reversed(['#FFAAAA', '#AAFFAA', '#AAAAFF'])))
    cmap_bold = ListedColormap(list(reversed(['#FF0000', '#00FF00', '#0000FF'])))
    
    # Plot the decision boundary. For that, we will assign a color to each
    # point in the mesh [x_min, m_max]x[y_min, y_max].
    h = .05  # step size in the mesh
    x_min, x_max = X.values[:, 0].min() - .5, X.values[:, 0].max() + .5
    y_min, y_max = X.values[:, 1].min() - .5, X.values[:, 1].max() + .5
    xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
    Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])

    # Put the result into a color plot
    Z = Z.reshape(xx.shape)
    
    if fig is None:
        figura = plt.figure(figsize=figsize)
        plt.pcolormesh(xx, yy, Z, cmap=cmap_light, alpha=0.4)

        plt.scatter(X.values[:, 0], X.values[:, 1], c=y.values,
                    cmap=cmap_bold, edgecolors='k')
        plt.title(title)
        plt.xlabel(xlabel)
        plt.ylabel(ylabel)
        #plt.xticks(())
        #plt.yticks(())
        plt.xlim(xx.min(), xx.max())
        plt.ylim(yy.min(), yy.max())
        
    else:
        fig.pcolormesh(xx, yy, Z, cmap=cmap_light, alpha=0.4)

        fig.scatter(X.values[:, 0], X.values[:, 1], c=y.values,
                    cmap=cmap_bold, edgecolors='k')
        fig.set_title(title)
        fig.set_xlabel(xlabel)
        fig.set_ylabel(ylabel)
        fig.set_xticks(())
        fig.set_yticks(())
        fig.set_xlim(xx.min(), xx.max())
        fig.set_ylim(yy.min(), yy.max())

No es la función más bonita ni agradable de leer del mundo. Pero vamos a probarla para dibujar nuestra clasificación en los datos de training y testing:

In [28]:
figura = plt.figure(figsize=(16,7))
figura.suptitle("Regresión logística", fontsize=14, fontweight="bold") # Supertítulo, común a ambos subplots

# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(lr,
                      dataset_clasificacion_training[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_training["label"],
                      title="Training",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes1)

# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(lr,
                      dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_testing["label"],
                      title="Testing",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes2)

Así es como nuestra regresión logística ha aprendido a separar los casos de calvo frente a los de no calvo. Aprendió a hacerlo en el training set, y hemos medido su rendimiento en el testing set. Vemos que la regresión logística separa linealmente los datos (con una línea recta). Puede que si te fijes bien veas unos pequeños escaloncitos en el dibujo de la línea recta, pero es tema de la visualización: la línea de separación es recta en la realidad.

Podemos ver también que en el conjunto de testing no acierta a predecir todos los puntos (hay puntos rojos en la zona donde el modelo decide que esos son azules, y viceversa). Aún así, hemos comprobado con lr.score(features_testing, labels_testing) que este error $E_{out}$ es menor al 10%. El modelo ni siquiera es capaz de clasificar bien todos los puntos en el dataset de entrenamiento... ($E_{in}$ no es 0%).

No obstante, la regresión logística es uno de los clasificadores más utilizados. Es sencillo de entender y de implementar, y es bastante rápido entrenando.

Muy rápido, vamos a probar a tocar alguno de los (hiper)parámetros de la regresión logística, para ver si afecta mucho al resultado. Para ello, creamos una nueva:

In [29]:
lr2 = LogisticRegression(penalty="l2", # Penalización L2 o Ridge
                         solver="lbfgs", # El algoritmo que se usa para minimizar el error, y entrenar el modelo
                         max_iter=600) # Número máximo de iteraciones para el entrenamiento 

Entrenamos con el conjunto de training:

In [30]:
lr2.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
       y=dataset_clasificacion_training["label"])
Out[30]:
LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,
          intercept_scaling=1, max_iter=600, multi_class='ovr', n_jobs=1,
          penalty='l2', random_state=None, solver='lbfgs', tol=0.0001,
          verbose=0, warm_start=False)

Comprobamos el rendimiento en el conjunto de testing:

In [31]:
lr2.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
         y=dataset_clasificacion_testing["label"])
Out[31]:
0.94166666666666665

Y dibujamos:

In [32]:
figura = plt.figure(figsize=(16,7))
figura.suptitle("Regresión logística (regularización L2, L-BFGS, 600 iteraciones)", fontsize=14, fontweight="bold") # Supertítulo, común a ambos subplots

# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(lr2,
                      dataset_clasificacion_training[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_training["label"],
                      title="Training",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes1)

# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(lr2,
                      dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_testing["label"],
                      title="Testing",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes2)

No ha cambiado mucho la cosa (algo ha mejorado, sí). En problemas de 2 dimensiones (2 variables/features, en este caso hormona_a y hormona_b), las posibles regularizaciones que podamos elegir no suelen tener mucho efecto en la regresión logística (especialmente si los datos no tienen valores extremos o mucho ruido).

Por ello, vamos a seguir con el mismo dataset y probar otros clasificadores.

3.D.b. Árboles

Los árboles de decisión son muy sencillos de entender, rápidos de entrenar (cuando es de uno en uno; luego cuando veamos ensembles veremos que no son tan rápidos) y son modelos no lineales, es decir: que la línea que separe los datos no tiene que ser necesariamente recta.

Importemos el modelo y entrenémoslo. El modelo es el sklearn.tree.DecisionTreeClassifier:

In [33]:
from sklearn.tree import DecisionTreeClassifier

arbol1 = DecisionTreeClassifier(criterion="entropy", # Criterio de división (entropía o Gini soportado)
                                max_depth=5, # Máxima profundidad de árbol
                                min_samples_split=1, # Número mínimo de puntos en un nodo para hacer split
                                min_samples_leaf=1 # Número mínimo de puntos para considerar un nodo un nodo-hoja
                                )
In [34]:
arbol1.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
           y=dataset_clasificacion_training["label"])
Out[34]:
DecisionTreeClassifier(class_weight=None, criterion='entropy', max_depth=5,
            max_features=None, max_leaf_nodes=None,
            min_impurity_split=1e-07, min_samples_leaf=1,
            min_samples_split=1, min_weight_fraction_leaf=0.0,
            presort=False, random_state=None, splitter='best')

Entrenado. Vamos a ver qué rendimiento nos da en el conjunto de testing:

In [35]:
arbol1.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
             y=dataset_clasificacion_testing["label"])
Out[35]:
0.90833333333333333

Y dibujemos:

In [36]:
figura = plt.figure(figsize=(16,7))
figura.suptitle("Arbol de decisión (prof. max. 5)", fontsize=14, fontweight="bold") # Supertítulo, común a ambos subplots

# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(arbol1,
                      dataset_clasificacion_training[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_training["label"],
                      title="Training",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes1)

# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(arbol1,
                      dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_testing["label"],
                      title="Testing",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes2)

Ahora vemos que las zona de decisión es no lineal (esa forma de "rectangulitos" es muy peculiar de los árboles). No obstante, nuestro rendimiento en el conjunto de testing es peor que en la regresión logística, a pesar de que el modelo es más complejo. ¿Por qué?

Porque estamos sobreajustando: el árbol "se está aprendiendo demasiado" los datos de entrenamiento, y cuando tiene que clasificar los del conjunto de testing, comete más fallos que si separara los datos con una sencilla línea recta. Es el principal inconveniente de los árboles: tienden a sobreajustar bastante. Cuando veamos los ensembles veremos distintas formas de evitar que los árboles sobreajusten.

A la gente le suele gustar "ver" los árboles, para ver las reglas de decisión. Sklearn nos permite visualizarlo fácilmente, utilizando sklearn.tree.export_graphviz. La visualización se exporta en formato DOT, el cual podemos pasar a pdf, imagen png y otros. Vamos a ver cómo pasarlo a png y mostrarlo en el notebook.

Nota: antes de instalar las bibliotecas mostradas en la siguiente celda, es posible que necesites instalar graphviz en tu sistema operativo.

In [37]:
from sklearn.tree import export_graphviz

# También necesitamos la biblioteca pydot si usamos Python 2,
# o pydotplus si usamos Python 3. Estas bibliotecas permiten gestionar
# los archivos de formato DOT. Puedes instalar pydotplus con pip:
import pydotplus

# Además, vamos a usar un truco: en vez de exportar el archivo dot
# como archivo al disco duro, vamos a meterlo en una variable
# utilizando la siguiente biblioteca para escribir el output
# del archivo DOT a un string:
from sklearn.externals.six import StringIO

# Creamos una instancia de StringIO, que es un string
# especial donde vamos a "volcar" el archivo DOT:
archivo_dot_en_string = StringIO()

# Exportamos el gráfico:
export_graphviz(arbol1,
                out_file=archivo_dot_en_string,
                feature_names=["hormona_a", "hormona_b"],
                class_names=["no_calvo","calvo"], # Por orden de menor a mayor. 0 es no_calvo, y 1 es calvo
                filled=True, # Con colores
                rounded=True # Con esquinas redondeadas
               )

# Usamos pydotplus para procesar el archivo DOT:
grafo_del_arbol = pydotplus.graph_from_dot_data(archivo_dot_en_string.getvalue())

# Y ahora lo podemos mostrar en el notebook con IPython.display.Image:
from IPython.display import Image
Image(grafo_del_arbol.create_png(), width=900)
Out[37]:

Mucho split para solo dos features... Apesta a sobreajuste. Vamos a lanzar uno con menos profundidad:

In [38]:
arbol2 = DecisionTreeClassifier(criterion="entropy",
                                max_depth=3,
                                min_samples_split=5,
                                min_samples_leaf=5
                                )
arbol2.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
           y=dataset_clasificacion_training["label"])

# Score en testing:
arbol2.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
             y=dataset_clasificacion_testing["label"])
Out[38]:
0.92500000000000004

Ahora tiene menor $E_{out}$.

In [39]:
figura = plt.figure(figsize=(16,7))
figura.suptitle("Arbol de decisión (prof. max. 3)", fontsize=14, fontweight="bold") # Supertítulo, común a ambos subplots

# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(arbol2,
                      dataset_clasificacion_training[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_training["label"],
                      title="Training",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes1)

# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(arbol2,
                      dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_testing["label"],
                      title="Testing",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes2)
In [40]:
archivo_dot_en_string = StringIO()
export_graphviz(arbol2,
                out_file=archivo_dot_en_string,
                feature_names=["hormona_a", "hormona_b"],
                class_names=["no_calvo","calvo"],
                filled=True,
                rounded=True
               )
grafo_del_arbol = pydotplus.graph_from_dot_data(archivo_dot_en_string.getvalue())
Image(grafo_del_arbol.create_png(), width=900)
Out[40]:

Mucho mejor. Pasemos al siguiente clasificador.

3.D.c. Nearest Neighbors

Nearest Neighbors (o K-Nearest Neighbors o KNN) es también un modelo muy sencillo de entender. También es no lineal (es decir, que genera regiones de decisión no lineales), y en algunos casos da un gran balance entre sencillez de modelo y resultados. Se basa en un principio muy básico: si conocemos el label de un determinado punto, los puntos cercanos deberían tener el mismo label; y esa será su predicción. Dependiendo de cuántos vecinos contemplemos para un punto a clasificar, el resultado será distinto.

Por ejemplo: supongamos 5-NN, o contemplar los 5 puntos conocidos más cercanos al que queremos clasificar. Si tenemos ese punto a clasificar rodeado de tres calvos y dos no-calvos, se predicirá que el punto a clasificar es calvo (porque de sus 5 vecinos más cercanos, la mayoría (tres) son calvos).

Ahora, para el mismo punto a clasificar, contemplamos 7-NN; y resulta que tres conocidos son calvos (los mismos 3 de antes), pero los otros cuatro conocidos más cercanos son no-calvos (los dos de antes, más otros dos). Ahora, el punto será predicho como no-calvo.

Nearest Neighbors es un modelo no lineal, como veremos en los plots de más abajo.

Sklearn implementa el clasificador KNN a través de sklearn.neighbors.KNeighborsClassifier. Probémoslo:

In [41]:
from sklearn.neighbors import KNeighborsClassifier

knn1 = KNeighborsClassifier(n_neighbors=5, # Número de vecinos contemplados
                            metric="euclidean", # Distancia utilizada para ver cuáles son los más cercanos
                            algorithm="brute" # Algoritmo de búsqueda de vecinos (por fuerza bruta)
                           )

knn1.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
           y=dataset_clasificacion_training["label"])
Out[41]:
KNeighborsClassifier(algorithm='brute', leaf_size=30, metric='euclidean',
           metric_params=None, n_jobs=1, n_neighbors=5, p=2,
           weights='uniform')

Rendimiento en el conjunto testing:

In [42]:
knn1.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
             y=dataset_clasificacion_testing["label"])
Out[42]:
0.94166666666666665

Plots:

In [43]:
figura = plt.figure(figsize=(16,7))
figura.suptitle("5-NN", fontsize=14, fontweight="bold") # Supertítulo, común a ambos subplots

# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(knn1,
                      dataset_clasificacion_training[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_training["label"],
                      title="Training",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes1)

# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(knn1,
                      dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_testing["label"],
                      title="Testing",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes2)

No está mal. Probemos con 1-NN (es decir, contemplar solo el punto más cercano):

In [44]:
knn2 = KNeighborsClassifier(n_neighbors=1, # Número de vecinos contemplados
                            metric="euclidean", # Distancia utilizada para ver cuáles son los más cercanos
                            algorithm="brute" # Algoritmo de búsqueda de vecinos (por fuerza bruta)
                           )

knn2.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
           y=dataset_clasificacion_training["label"])

print("Bien clasificados en testing: " + str(knn2.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
                                 y=dataset_clasificacion_testing["label"])))

# Plot:
figura = plt.figure(figsize=(16,7))
figura.suptitle("1-NN", fontsize=14, fontweight="bold") # Supertítulo, común a ambos subplots

# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(knn2,
                      dataset_clasificacion_training[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_training["label"],
                      title="Training",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes1)

# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(knn2,
                      dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_testing["label"],
                      title="Testing",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes2)
Bien clasificados en testing: 0.891666666667

Mayor $E_{out}$, y sobreajuste obvio.

El resultado de K-NN puede cambiar ampliamente en función de la métrica de distancia que utilicemos. La euclídea funciona bien para espacios dimensionales pequeños (este por ejemplo es 2-D); pero cuando tenemos muchas features, es mejor utilizar otra de las métricas de distancia que nos ofrece sklearn.

3.D.d. Naïve Bayes

Naïve Bayes son una familia de modelos puramente probabilísticos que son bastante utilizados en clasificación. Como el propio nombre indica, están fundamentados en el teorema de Bayes. Dada una label (en nuestro caso calvo/no-calvo) y una serie de features (en nuestro caso hormona_a y hormona_b), Bayes dice que la relación es:

$$P(calvo \mid hormona_a, hormona_b) = \frac{P(calvo) P(hormona_a, hormona_b \mid calvo)}{P(hormona_a, hormona_b)}$$

Dado esto, Naïve Bayes asume que no existe ninguna relación entre la concentración de la hormona a y la de la hormona b. Esto es, que son probabilísticamente independientes (su probabilidad condicionada es la misma que sus probabilidades independientes). Esto es lo que hace Naïve (o ingenuo) al modelo: en la realidad, sí suelen existir relaciones entre las features; sean conocidas o no.

Dentro de sklearn.naive_bayes tenemos varios modelos, en función de cómo asumimos que es la verosimilitud de las features o evidencias, es decir, cómo creemos que es $P(x_i \mid y)$.

Si bien muchos de los otros modelos de machine learning están basados en mayor o menor medida en verosimilitud y modelos probabilísticos, Naïve Bayes está puramente basado en ello. El principal problema de los modelos probabilísticos "puros" es que nos obligan a hacer asunciones sobre cuál es la distribución de probabilidad de las features; y todo lo basado en asunciones, puede ser muy erróneo...

Probemos pues. Vamos a usar un GaussianNB:

In [45]:
from sklearn.naive_bayes import GaussianNB

nb = GaussianNB() # GaussianNB pasa del tema argumentos y parámetros de regularización...

nb.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
           y=dataset_clasificacion_training["label"])

print("Bien clasificados en testing: %.15f" % nb.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
                                 y=dataset_clasificacion_testing["label"]))

# Plot:
figura = plt.figure(figsize=(16,7))
figura.suptitle("Naïve Bayes", fontsize=14, fontweight="bold") # Supertítulo, común a ambos subplots

# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(nb,
                      dataset_clasificacion_training[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_training["label"],
                      title="Training",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes1)

# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(nb,
                      dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_testing["label"],
                      title="Testing",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes2)
Bien clasificados en testing: 0.950000000000000

Dependiendo de la verosimilitud de las features (en este caso gausiana, de ahí lo de GaussianNB), Naïve Bayes será lineal o no lineal. Por lo general, nos encontraremos que el Naïve Bayes es un clasificador no lineal (como en el ejemplo que acabamos de hacer, aunque no se aprecia demasiado).

3.D.e. Otros clasificadores: SVMs, redes y demás

Sklearn incluye otros muchos clasificadores, y sobre todo, quasi-infinitas variaciones de los ya vistos (gracias a la cantidad de hiperparámetros que podemos tocar para cada modelo). Por último, vamos a probar los SVMs (Support Vector Machines), los cuales se verán en profundidad en el módulo de algoritmos avanzados.

Siendo perfectamente discutible: desde mi experiencia, los SVMs son los mejores clasificadores que hay (más que las redes neuronales, ensembles de árboles, y demás), y la matemática detrás de su formulación es una belleza.

Probemos rápidamente un SVM:

In [46]:
from sklearn.svm import SVC

svm = SVC(kernel="rbf") # Kernel de función base radial (ya lo verás)

svm.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
           y=dataset_clasificacion_training["label"])

print("Bien clasificados en testing: %.16f" % svm.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
                                 y=dataset_clasificacion_testing["label"]))


# Plot:
figura = plt.figure(figsize=(16,7))
figura.suptitle("SVM con kernel RBF", fontsize=14, fontweight="bold") # Supertítulo, común a ambos subplots

# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(svm,
                      dataset_clasificacion_training[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_training["label"],
                      title="Training",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes1)

# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(svm,
                      dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_testing["label"],
                      title="Testing",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes2)
Bien clasificados en testing: 0.9500000000000000

Desde la versión 0.18 (28 de septiembre de 2016) Sklearn incluye tres tipos de redes:

  • La RBM (Restricted Boltzmann machines), que no sirven para clasificación como tal, en su formulación pura (se utilizan sobre todo como autoencoders). Las RBMs son muy buenas aprendiendo distribuciones de probabilidad. Este post de los señores de DeepLearning4J explica muy bien cómo funciona una RBM. La RBM ya estaba disponible en versiones anteriores de sklearn. Las otras dos son nuevas de la versión 0.18.

  • MLPClassifier implementa la versión clasificador del MLP (Perceptrón Multicapa). El perceptrón multicapa es la red neuronal más sencilla que hay, pero aún así es la más utilizada para datasets convencionales (los que no son de imágenes, texto, secuencias, vídeos... Sino variables "de negocio"). Podemos entender el perceptrón multicapa como una sucesión de capas de regresiones logísticas (más o menos).

  • MLPRegressor implementa el MLP pero como regresor en vez de clasificador (de forma que la última capa es como las regresiones lineales).

En general, el diseño de redes neuronales requiere tocar muchos más hiperparámetros que el resto de modelos que hemos visto. De ahí que existan bibliotecas específicas para el diseño y entrenamiento de redes neuronales... Y la mayoría de ellas están escritas para Python. Destaco:

  • Theano: la de los campeones. Llamarla solo como una biblioteca para hacer redes es casi insultar a Theano. Es una biblioteca capaz de manejar tensores (piensa en arrays de Numpy) de forma simbólica, para delegar la ejecución en tarjetas gráficas o GPUs, que son mucho más rápidas para entrenar redes. Muy manual y complicada: una red sencilla son unas 1.500 líneas de código fácilmente.
  • TensorFlow: hecha principalmente por Google, es más sencilla de utilizar. No es mi favorita en absoluto (soy más the Theano), pero va en gustos. Cuando se liberó hace aproximadamente un año era muy pobre, tanto en rendimiento como en funcionalidad. Pero va mejorando, y actualmente existe una biblioteca (también parte del proyecto) llamada TFLearn que es similar a Keras. Independientemente de que la acabes usando o no, existe mucha documentación sobre ella asociada a aprender sobre redes.
  • Keras: probablemente la más recomendable para empezar (siempre y cuando ya sepas de redes, y solo te quede por conocer una buena bilioteca para empezar a entrenarlas). Keras es una capa de abstracción por encima tanto de Theano como de TensorFlow. Esto quiere decir: diseñas redes de forma relativamente sencilla con Keras, y se ejecutan utilizando los motores bien the Theano, bien de TensorFlow (recomiendo encarecidamente utilizar el de Theano a día de hoy). Está teniendo un crecimiento brutal en los últimos meses, y es la favorita de muchos.

Ejemplo del MLP de sklearn 0.18 en acción (me tocó usarlo en el curro para una ppt...). Es curioso ver cómo va mejorando el ajuste con las iteraciones: MLP

Vamos a entrenar también un MLP con nuestro dataset (aunque el efecto es menos vistoso, básicamente porque nuestro dataset se puede separar casi con una línea recta):

In [47]:
from sklearn.neural_network import MLPClassifier

# Los hiperparámetros de las redes son difíciles de entender sin haber estudiado bien las redes.
# En el módulo de algoritmos avanzados veréis cómo funcionan estos parámetros:

mlp = MLPClassifier(hidden_layer_sizes=(15,15), # Número de neuronas en las capas ocultas
                    activation="logistic", # En cada neurona se aplica la función logística
                    solver="sgd", # Descenso gradiente estocástico
                    batch_size=10, # Cuán "estocástico" es el descenso gradiente (10 observaciones por "epoch")
                    alpha=0.0,
                    learning_rate="adaptive",
                    learning_rate_init=0.015,
                    momentum=0.8,
                    nesterovs_momentum=True,
                    max_iter=500)

mlp.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
           y=dataset_clasificacion_training["label"])

print("Bien clasificados en testing: %.16f" % mlp.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
                                                        y=dataset_clasificacion_testing["label"]))


# Plot:
figura = plt.figure(figsize=(16,7))
figura.suptitle("MLP (dos capas ocultas de 15 neuronas),\nactivación logística/sigmoide", fontsize=14, fontweight="bold") # Supertítulo, común a ambos subplots

# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(mlp,
                      dataset_clasificacion_training[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_training["label"],
                      title="Training",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes1)

# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(mlp,
                      dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_testing["label"],
                      title="Testing",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes2)
Bien clasificados en testing: 0.9416666666666667

Ahora continuaremos con los clasificadores, pero de otra manera. Vamos a empezar con los ensembles.

4. Ensemble Learning

El principio del Ensemble Learning es sencillo: la unión hace la fuerza. Consiste en combinar varios modelos en uno solo, con el fin de que el modelo final sea más robusto y generalice mejor.

La intuición es la siguiente: mediante la aplicación de varios modelos a la vez, los fallos de uno de los modelos son compensados por los aciertos de otros de los modelos sobre los mismos puntos. Los ensembles son muy utilizados en la práctica, y la mayoría de los Kaggle se ganan mediante la generación de ensembles de modelos.

Sklearn implementa varias técnicas de ensembles, que se dividen básicamente en dos grupos:

  • Bagging, donde se lanzan varios modelos a la vez (normalmente cada uno de ellos con subconjunto del training set). Las predicciones finales se basan en el consenso al que llegan los modelos lanzados: se puede poner a los modelos a votar en el caso de clasificadores, hacer la media de la predicción en el caso de regresores...
  • Boosting, donde primero se lanza un modelo. Luego, se lanza otro utilizando parte de las observaciones bien clasificadas por el primero, y las observaciones mal clasificadas por éste (o modificaciones sobre dichos datos). Y así sucesivamente. Se podría decir que el segundo modelo intenta corregir los fallos del primero, el tercero intenta corregir los del segundo, etcétera.

sklearn.ensemble es el lugar donde se encuentran los ensembles. Algunos de ellos son los siguientes:

Técnica Módulo Comentario
Bagging ensemble.BaggingClassifier Media/votación de múltiples ejecuciones del estimador sobre diferentes subconjuntos de datos
Random Forests ensemble.RandomForestsClassifier Bagging de árboles de decisión
Extremelly Randomized Trees ensemble.ExtraTreesClassifier Lanza árboles randomizados
ADABoost ensemble.AdaBoostClassifier Adaptative Boosting
Gradient Boosted Trees ensemble.GradientBoostingClassifier Boosting de árboles

Vamos a explorar algunos:

4.A. Bagging

Vamos a probar a hacer un bagging de regresiones logísticas sobre nuestro dataset sobre la alopecia, utilizando el BaggingClassifier:

In [48]:
from sklearn.ensemble import BaggingClassifier

# Generemos nuestra regresión logística:
lr = LogisticRegression(penalty="l1", # Penalización L1 o LASSO
                        solver="liblinear")

# Y ahora generamos el bagging:
bagging_lr = BaggingClassifier(lr, # Modelos
                               n_estimators=50, # Número de modelos en el bagging
                               max_samples=0.3, # Proporción máxima del dataset dada a cada modelo para entrenar
                               max_features=1,  # Número máximo de features dado a cada modelo para entrenar
                               bootstrap=True,  # Con reemplazamiento en la parte del dataset cogida por modelo
                               bootstrap_features=True, # Con reemplazamiento en las features
                               n_jobs=4         # Paralelizar el trabajo en 4 núcleos
                              )

Al igual que con un modelo normal y corriente, lo entrenamos con .fit():

In [49]:
bagging_lr.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
               y=dataset_clasificacion_training["label"])
Out[49]:
BaggingClassifier(base_estimator=LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,
          intercept_scaling=1, max_iter=100, multi_class='ovr', n_jobs=1,
          penalty='l1', random_state=None, solver='liblinear', tol=0.0001,
          verbose=0, warm_start=False),
         bootstrap=True, bootstrap_features=True, max_features=1,
         max_samples=0.3, n_estimators=50, n_jobs=4, oob_score=False,
         random_state=None, verbose=0, warm_start=False)

Y comprobamos el rendimiento en el conjunto de testing:

In [50]:
print("Bien clasificados en testing: %.16f" % bagging_lr.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
                                                               y=dataset_clasificacion_testing["label"]))
Bien clasificados en testing: 0.9500000000000000
In [51]:
figura = plt.figure(figsize=(16,7))
figura.suptitle("Bagging de LR (50 modelos)", fontsize=14, fontweight="bold")

# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(bagging_lr,
                      dataset_clasificacion_training[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_training["label"],
                      title="Training",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes1)

# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(bagging_lr,
                      dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_testing["label"],
                      title="Testing",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes2)

Vemos que lo hace muy bien en el conjunto de testing. Aunque no se aprecia, el resultado es no lineal. Parece que el modelo generaliza muy bien.

4.B. Random Forests

Random Forests es uno de los modelos más famosos para clasificación. No deja de ser un bagging de árboles de decisión, con una pequeña modificación (básicamente que cada árbol es entrenado con subconjuntos no solo de las observaciones, sino también de las features). Pero un bagging de árboles al fin y al cabo. Suele dar muy buenos resultados.

Probemos entonces el RandomForestClassifier:

In [52]:
from sklearn.ensemble import RandomForestClassifier

rf = RandomForestClassifier(n_estimators=20, # 20 árboles
                            criterion="entropy", # Igual que en el DecisionTreeClassifier
                            max_depth=3, # Igual que en el DecisionTreeClassifier
                            min_samples_split=10,
                            min_samples_leaf=5,
                            bootstrap=True)

rf.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
       y=dataset_clasificacion_training["label"])

print("Bien clasificados en testing: %.16f" % rf.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
                                                       y=dataset_clasificacion_testing["label"]))
Bien clasificados en testing: 0.9333333333333333
In [53]:
figura = plt.figure(figsize=(16,7))
figura.suptitle("Random Forests (20 árboles)", fontsize=14, fontweight="bold")

# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(rf,
                      dataset_clasificacion_training[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_training["label"],
                      title="Training",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes1)

# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(rf,
                      dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_testing["label"],
                      title="Testing",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes2)

4.C. AdaBoost

AdaBoost (o Adaptive Boosting) es el algoritmo de boosting más popular y, en muchas ocasiones, el que mejores resultados da.

AdaBoost va pasando clasificadores (o regresores) sobre formas modificadas de los datos, de forma repetida. Por cada boosting iteration (pasada de un nuevo clasificador), se asignan pesos o importancias a cada punto, de forma que los puntos mal clasificados por el modelo anterior cobran más importancia; mientras que los puntos bien clasificados la van perdiendo. Esta técnica adaptativa de boosting es muy popular y bastante efectiva en general.

Utilicemos sklearn.ensemble.AdaBoostClassifier para generar un adaptive boosting de regresiones logísticas:

In [54]:
from sklearn.ensemble import AdaBoostClassifier

lr = LogisticRegression(penalty="l2",
                        solver="lbfgs")

ab = AdaBoostClassifier(base_estimator=lr,
                        n_estimators=20, # Número de clasificadores
                        learning_rate=0.4 # La "contribución" de cada clasificador 
                        )

ab.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
       y=dataset_clasificacion_training["label"])

print("Bien clasificados en testing: %.16f" % ab.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
                                                       y=dataset_clasificacion_testing["label"]))
Bien clasificados en testing: 0.9583333333333334
In [55]:
figura = plt.figure(figsize=(16,7))
figura.suptitle("AdaBoost (Logistic Regression)", fontsize=14, fontweight="bold")

# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(ab,
                      dataset_clasificacion_training[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_training["label"],
                      title="Training",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes1)

# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(ab,
                      dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_testing["label"],
                      title="Testing",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes2)

4.D. Gradient Boosting

La técnica de Gradient Boosting es, cuanto menos, curiosa: al igual que el resto de técnicas de boosting, se lanza primero un clasificador (en el caso de Gradient Boosting suelen ser árboles de decisión). Luego, se formula el árbol generado como una función diferenciable; y optimizable (mediante la computación de gradientes e ir descendiendo el error generado) a través del árbol generado para la siguiente iteración. Los señores de XGBoost (biblioteca que mencionaremos al final) lo explican mejor en su blog.

De hecho, XGBoost es una implementación de Gradient Boosting (no incluida en scikit-learn, sino como una biblioteca de Python aparte) escalable y que incorpora una buena regularización (luego hablaremos de las regularizaciones). Pero aquí vamos a utilizar la implementación de Gradient Boosting de scikit-learn.

In [56]:
from sklearn.ensemble import GradientBoostingClassifier

gbt = GradientBoostingClassifier(learning_rate=0.05,
                                 n_estimators=40,
                                 max_depth=3, # Son árboles al fin y al cabo
                                 min_samples_split=5,
                                 min_samples_leaf=5
                                 )

gbt.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
       y=dataset_clasificacion_training["label"])

print("Bien clasificados en testing: %.16f" % gbt.score(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
                                                       y=dataset_clasificacion_testing["label"]))
Bien clasificados en testing: 0.9333333333333333
In [57]:
figura = plt.figure(figsize=(16,7))
figura.suptitle("Gradient Boosted Trees (30 Trees)", fontsize=14, fontweight="bold")

# Dibujamos en el dataset de training:
axes1=plt.subplot(1,2,1)
plot_decision_regions(gbt,
                      dataset_clasificacion_training[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_training["label"],
                      title="Training",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes1)

# Dibujamos en el dataset de testing:
axes2=plt.subplot(1,2,2)
plot_decision_regions(gbt,
                      dataset_clasificacion_testing[["hormona_a", "hormona_b"]],
                      dataset_clasificacion_testing["label"],
                      title="Testing",
                      xlabel="Concentración de hormona a",
                      ylabel="Concentración de hormona b",
                      fig=axes2)

Los ensembles de modelos se utilizan mucho en clasificación; pero también se pueden utilizar en problemas de regresión: apartado que veremos después de las métricas de clasificación.

5. Métricas de clasificación

Hasta ahora hemos utilizado una métrica para medir el rendimiento de nuestros clasificadores. Al aplicar el método .score() en un clasificador entrenado, nos devuelve un número entre 0 y 1. Este número nos dice el porcentaje de observaciones (o puntos) bien clasificadas; métrica conocida como accuracy.

Esta métrica está muy bien en algunos casos, pero no siempre. Supongamos un dataset como el siguiente (vamos a crear solo las labels, por simplificar):

In [58]:
np.concatenate((np.zeros(90),np.ones(10)))
Out[58]:
array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,
        1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.])

Donde tenemos 90 observaciones que son de la clase cero, y solo 10 que son de la clase 1. Pues bien, supongamos que creamos un "modelo" como el siguiente:

Nuestro modelo predice siempre 0, independientemente de la observación.

Con este modelo tan estúpido conseguiríamos un 90% de puntos correctamente clasificados. No obstante, eso no quiere decir que sea un buen modelo (evidentemente no lo es). Más importante que clasificar bien un X% de los puntos es que nuestro modelo sea capaz de discriminar, y diferenciar lo que debería ser un uno de lo que debería ser un cero.

Por esa razón la métrica de porcentaje de puntos bien clasificados no es siempre ideal.

5.A. Confusion matrix

Afortunadamente, existen otras métricas que nos permiten ver el rendimiento de un clasificador.

Una de las primeras cosas que podemos hacer tras entrenar un clasificador es construir lo que se conoce como confusion matrix o matriz de confusión. La matriz de confusión no es más que una tablita que se utiliza normalmente para ver qué tal lo ha hecho un clasificador en problemas de clasificación binaria, es decir, problemas de clasificar dos labels (aunque se puede extender a problemas de clasificar más de dos clases). Esta tablita tiene cuatro cuadrantes en el caso de la clasificación binaria:

Predicho 1 Predicho 0
Valor real 1 True Positives False Negatives
Valor real 0 False Positives True Negatives
  • El cuadrante superior izquierdo muestra los true positives, es decir, el número de unos bien clasificados.
  • El cuadrante inferior izquierdo muestra los false positives: el número de ceros clasificados como uno.
  • El cuadrante superior derecho muestra los false negatives: el número de unos clasificados como cero.
  • El cuadrante inferior derecho muestra los true negatives: el número de ceros bien clasificados.

Esta matriz nos permite, de un rápido vistazo, determinar si estamos "sobreclasificando" de alguna de las dos clases (por ejemplo: muchísimos ceros clasificados como unos, pero pocos unos clasificados como ceros).

sklearn.metrics incluye la función confusion_matrix que nos permite construir la matriz de confusión. Toma como primer argumento las labels reales que vemos en las observaciones, y como segundo argumento las labels predichas por un modelo (las cuales podemos conseguir con modelo.predict()).

Dicho esto, vamos a construir la matriz de confusión del último clasificador que hemos entrenado, el gradient boosted trees que hemos llamado gbt. Utilizaremos los datos del conjunto de testing para ello:

In [59]:
from sklearn.metrics import confusion_matrix

matriz_confusion = confusion_matrix(y_true=dataset_clasificacion_testing["label"],
                                    y_pred=gbt.predict(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]])
                                   )

# Sklearn nos saca en la primera columna los
# predichos como cero, en vez de los predichos
# como uno, así que cambiamos el orden de las
# columnas y filas para que sea como el de nuestra
# definición previa:
matriz_confusion = np.array([(matriz_confusion[1]),
                              matriz_confusion[0]])[:,::-1]

matriz_confusion
Out[59]:
array([[55,  6],
       [ 2, 57]])

Acorde a los cuatro cuadrantes definidos antes, nuestro gbt en el conjunto de testing:

  • Ha predicho bien 55 unos -> true positives
  • Ha predicho 2 ceros como uno (2 hombres no calvos como calvos) -> false positives
  • Ha predicho 6 unos como cero (6 hombres calvos como no calvos) -> false negatives
  • Ha predicho bien 57 ceros -> true negatives

En el gráfico que hicimos cuando entrenamos el modelo resulta fácil de verificar.

Podemos ver que no existe ningún "fallo gordo" con el modelo. El problema vendría por ejemplo si las cantidades de false positives y de false negatives fueran muy dispares.

5.B. Precision, recall y F1-score

La matriz de confusión es fácil de interpretar, pero para nosotros; no para los ordenadores. Ellos prefieren un número, o un par de ellos más fáciles de comparar entre modelos.

Para eso existen precisamente otras métricas como son la precision, el recall y el F1-score:

5.B.a. Precision

La precision puede entenderse como la habilidad que tiene un clasificador para no predecir los ceros como unos. Su cálculo es muy sencillo:

$$precision=\frac{\sum true\,positives}{\sum (true\,positives + false\,positives)}$$

O dicho de otro modo:

$$precision=\frac{\sum Los\,adecuadamente\,predichos\,como\,1}{\sum Todos\,los\,predichos\,como\,1}$$

Podemos computar fácilmente la precision con sklearn.metrics.precision_score. Se utiliza exactamente igual que confusion_matrix, pero en este caso devuelve un solo número en vez de una matriz:

In [60]:
from sklearn.metrics import precision_score

precision_gbt = precision_score(y_true=dataset_clasificacion_testing["label"],
                                y_pred=gbt.predict(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]])
                               )

print("Precision: " + str(precision_gbt))
print("Cuando nuestro clasificador predice una observación como uno, acierta en el "+ str(precision_gbt*100)
      + "% de los casos (en el conjunto de testing).")
Precision: 0.964912280702
Cuando nuestro clasificador predice una observación como uno, acierta en el 96.4912280702% de los casos (en el conjunto de testing).

5.B.b. Recall

El recall es la habilidad que tiene un clasificador para predecir todos los unos reales. Su cálculo es:

$$recall=\frac{\sum true\,positives}{\sum (true\,positives + false\,negatives)}$$

O bien:

$$recall=\frac{\sum Los\,adecuadamente\,predichos\,como\,1}{\sum Todos\,los\,que\,son\,1\,en\,realidad}$$

Computamos el recall con sklearn.metrics.recall_score:

In [61]:
from sklearn.metrics import recall_score

recall_gbt = recall_score(y_true=dataset_clasificacion_testing["label"],
                          y_pred=gbt.predict(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]])
                         )

print("Recall: " +  str(recall_gbt))
print("Nuestro clasificador ha conseguido predecir bien el " + str(recall_gbt*100) 
      + "% de los unos (en el cojunto de testing).")
Recall: 0.901639344262
Nuestro clasificador ha conseguido predecir bien el 90.1639344262% de los unos (en el cojunto de testing).

5.B.c. F1-score

El F1-score esla métrica más utilizada dentro de las que son llamadas F-measures (valores-F en castellano). Básicamente combina en una sola métrica la precision y el recall a través de la siguiente fórmula:

$$F_1\mbox{-}score=2 \times\frac{precision \times recall}{precision + recall}$$

El F1-score da la misma importancia a la precision y al recall. Otras F-measures pueden utilizarse para dar más importancia a una métrica que a la otra. Qué F-measure utilizar depende mucho del negocio y del problema a resolver.

Utilicemos sklearn.metrics.f1_score para computarlo:

In [62]:
from sklearn.metrics import f1_score

f1_score_gbt = f1_score(y_true=dataset_clasificacion_testing["label"],
                        y_pred=gbt.predict(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]])
                       )

print("F1-score: " +  str(f1_score_gbt))
F1-score: 0.932203389831

Como regla general: cuanto más cercano sea el F1-score a 1, mejor; porque querrá decir que la precision y el recall también son cercanos a 1.

5.C. Curva ROC

La curva ROC (Receiver Operating Characteristic) es un gráfico muy útil para visualizar el rendimiento de un clasificador binario. Este gráfico tiene los siguientes ejes:

  • El eje x es el False Positive Rate (FPR).
  • El eje y es el True Positive Rate (TPR).

Donde el False Positive Rate es:

$$False\,Positive\,Rate=\frac{\sum false\,positives}{\sum (false\,positives+true\,negatives)}$$

O:

$$False\,Positive\,Rate=\frac{\sum Los\,0\,mal\,predichos\,como\,1}{\sum Todos\,los\,que\,son\,0\,en\,realidad}$$

Y el True Positive Rate es lo mismo que el recall:

$$True\,Positive\,Rate\,o\,recall=\frac{\sum true\,positives}{\sum (true\,positives + false\,negatives)}$$

Esta curva tiene dos usos principales:

  1. En el caso de clasificadores que generan un scoring, y que luego determinamos que son ceros las observaciones que puntúan por debajo de cierto valor umbral; y unos los que puntúan por encima: la curva ROC nos permite visualizar cuál es el mejor umbral que podemos poner. Por ejemplo: la regresión logística nos devuelve valores entre 0 y 1 para cada observación, y por defecto determinamos que los que den < 0.5 serán predichos como cero, y que los que den >= 0.5 serán predichos como uno. Pues bien, la curva ROC nos permite ver si existe otro umbral (0.7, 0.4...) que sea mejor para nuestro problema, en función de lo que busquemos. En estos casos, la curva ROC estará formada por los puntos de (FPR, TPR) para distintos umbrales.
  2. En el caso de clasificadores que no generen scoring la "curva ROC" será solo una línea que une el (0,0), el punto (FPR, TPR) de nuestro modelo, y el (1,1). Medir el área bajo esta curva (llamada Área Under ROC o AUC) es una métrica de rendimiento que nos permite comparar clasificadores (también se puede usar así para clasificadores que sí generan scoring).

Aquí nos centraremos en el segundo uso. Puedes aprender más sobre el primero con este vídeo.

Lo ideal es que un modelo tenga un False Positive Rate muy bajo (lo más cercano a $0$ posible), y un True Positive Rate muy elevado (lo más cercano a $1$ posible). De esta forma, el área bajo la curva ROC (AUC) será lo más elevada posible. Veámoslo con unos ejemplos.

Supongamos que hacemos un clasificador, y obtenemos que su FPR es $0.3$; y su TPR es $0.4$. Con este punto, el (0,0) y el (1,1) podemos dibujar la curva ROC:

In [63]:
coordenadas_fpr = [0, 0.3, 1]
coordenadas_tpr = [0, 0.4, 1]

plt.figure(figsize=(8,6))

# Primero pintamos el random guess:
plt.plot([0,1],[0,1], "g--")

# Y ahora el FPR y TPR de nuestro modelo:
plt.plot(coordenadas_fpr, coordenadas_tpr, "r")

# Coloreamos el área bajo la curva
plt.fill_between(coordenadas_fpr, coordenadas_tpr, alpha=0.4, color="yellow")

plt.title("Curva ROC (modelo mediocre)")
plt.xlabel("False Positive Rate")
plt.ylabel("True Positive Rate")
pass

La línea verde quebrada muestra lo que se llama como random guess o línea de no-discriminación, es decir: lo que haría un modelo que predice 0 o 1 "a boleo". Este modelo acertaría aproximadamente el 50% de los casos, dando por hecho que tenemos el mismo número de ceros que de unos. Se suele representar para comparar con la curva ROC de nuestro modelo.

Y la línea roja es la curva ROC de este modelo. Este modelo es poco mejor que uno que adivinara "a boleo": lo ideal es que la curva crezca lo más posible al principio, "despegándose" lo más posible de la línea de random guess por arriba. El área bajo esta curva ROC (coloreada de amarillo) es cercana a $0.5$. Teniendo en cuenta que el área total del gráfico es $1$, es poco más de la mitad: lo cual no es ideal. Fijándonos en el punto (0.3,0.4) de nuestro gráfico, podemos interpretar esta curva ROC como:

Para detectar y acertar el 40% de los unos, me "trago" o clasifico mal un 30% de los ceros de mi dataset.

Un clasificador perfecto tendría un FPR de $0$ y un TPR de $1$; y el área bajo su curva ROC sería $1$ (es decir, todo el gráfico).

Otro ejemplo: ahora supongamos que entrenamos otro modelo, y éste tiene un FPR también de 0.3; pero que por el contrario su TPR es más alto, siendo 0.8:

In [64]:
coordenadas_fpr = [0, 0.3, 1]
coordenadas_tpr = [0, 0.8, 1]

plt.figure(figsize=(8,6))

# Primero pintamos el random guess:
plt.plot([0,1],[0,1], "g--")

# Y ahora el FPR y TPR de nuestro modelo:
plt.plot(coordenadas_fpr, coordenadas_tpr, "r")

# Coloreamos el área bajo la curva
plt.fill_between(coordenadas_fpr, coordenadas_tpr, alpha=0.4, color="yellow")

plt.title("Curva ROC (modelo mejor)")
plt.xlabel("False Positive Rate")
plt.ylabel("True Positive Rate")
pass

Este está mejor: nuestra curva ROC ahora está bastante "más alta" que el random guess, y el AUC es bastante mejor. Sigue sin estar cerca de ser $1$, pero definitivamente es mayor; y bastante mayor que $0.5$.

Para detectar y acertar el 80% de los unos, me "trago" o clasifico mal un 30% de los ceros de mi dataset.

Y por último, supongamos un modelo muy malo, con FPR de 0.7 y TPR de 0.3:

In [65]:
coordenadas_fpr = [0, 0.7, 1]
coordenadas_tpr = [0, 0.3, 1]

plt.figure(figsize=(8,6))

# Primero pintamos el random guess:
plt.plot([0,1],[0,1], "g--")

# Y ahora el FPR y TPR de nuestro modelo:
plt.plot(coordenadas_fpr, coordenadas_tpr, "r")

# Coloreamos el área bajo la curva
plt.fill_between(coordenadas_fpr, coordenadas_tpr, alpha=0.4, color="yellow")

plt.title("Curva ROC (modelo penoso)")
plt.xlabel("False Positive Rate")
plt.ylabel("True Positive Rate")
pass

Este modelo es claramente pésimo: hace peor trabajo que el random guess, es decir: es mejor predecir "a boleo" que utilizar este modelo. Su AUC es bastante menor a $0.5$.

Para detectar y acertar el 30% de los unos, me "trago" o clasifico mal un 70% de los ceros de mi dataset.

Para usar un modelo como este, mejor nos dedicamos al arte abstracto.

Ahora, probemos a dibujar la curva ROC de nuestro modelo gbt.

sklearn.metrics.roc_curve nos genera las coordenadas para nuestros modelos de forma muy similar a como las hemos metido en los ejemplos que acabamos de inventarnos. La función nos devuelve una tupla con tres arrays de Numpy:

  1. El primer array nos da las coordenadas del eje x, es decir: del FPR.
  2. El segundo array nos da las coordenadas del eje y, es decir: el TPR.
  3. El tercer array nos da los thresholds o umbrales, en el caso de que quisiéramos utilizar la curva ROC para elegir el mejor umbral para un modelo de scoring (acorde al que pusimos como primero de los dos posibles usos para las curvas ROC). Así que aquí vamos a ignorarla: no tiene sentido utilizarla sin modelo de scoring.
In [66]:
from sklearn.metrics import roc_curve

roc = roc_curve(y_true=dataset_clasificacion_testing["label"],
                y_score=gbt.predict(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]]),
                pos_label=1.0 # Le decimos a roc_curve que los casos positivos son los unos.
               )

roc
Out[66]:
(array([ 0.        ,  0.03389831,  1.        ]),
 array([ 0.        ,  0.90163934,  1.        ]),
 array([2, 1, 0]))

Ahora cogemos las dos primeras arrays del resultado, y dibujamos la curva ROC de nuestro modelo gbt:

In [67]:
# Cogemos las coordenadas del resultado:
coordenadas_fpr = roc[0]
coordenadas_tpr = roc[1]

plt.figure(figsize=(8,6))

# Primero pintamos el random guess:
plt.plot([0,1],[0,1], "g--")

# Y ahora el FPR y TPR de nuestro modelo:
plt.plot(coordenadas_fpr, coordenadas_tpr, "r")

# Coloreamos el área bajo la curva
plt.fill_between(coordenadas_fpr, coordenadas_tpr, alpha=0.4, color="yellow")

plt.title("Curva ROC modelo gbt")
plt.xlabel("False Positive Rate")
plt.ylabel("True Positive Rate")
pass

Creo que queda claro que nuestro modelo gbt es muy bueno: la curva ROC se "despega" muchísimo del random guess por arriba; y crece muchísimo al principio. Por eso mismo, el AUC es casi toda la posible; y debe ser muy cercana a $1.0$.

De hecho, sklearn.metrics.roc_auc_score nos permite calcular fácilmente el AUC. Se utiliza como cualquier otra métrica de clasificación:

In [68]:
from sklearn.metrics import roc_auc_score

auc_gbt = roc_auc_score(y_true=dataset_clasificacion_testing["label"],
                        y_score=gbt.predict(X=dataset_clasificacion_testing[["hormona_a", "hormona_b"]])
                       )

print("AUC: " +  str(auc_gbt))
AUC: 0.933870519589

Efectivamente: el área bajo la curva ROC es muy cercana al máximo, que es $1$. Concluimos: AUC es una métrica sencilla de utilizar para evaluar la efectividad de nuestros clasificadores (al igual que la precision, recall y F1-score); y dibujar la curva ROC también es una forma vistosa de comparar clasificadores.

Puedes aprender más sobre la curva ROC aquí y aquí.

Existen muchas más métricas para medir el rendimiento de clasificadores, si bien las mencionadas son las más famosas y utilizadas. Otros ejemplos pueden ser la curva Lift (muy utilizada en marketing) o el LogLoss (disponible en sklearn).

Y hay todavía más, disponibles en sklearn.metrics: tanto para clasificación como para regresión y otros.

Ya está bien de clasificación. Pasemos al otro tipo principal de problema supervisado: la regresión.

6. Regresión

A diferencia de la clasificación (donde predecimos labels o clases), en regresión predecimos un número; también en función de una serie de features de cada observación. Un ejemplo podría ser el visto al principio de este documento: intentar predecir el precio de las casas en Boston a partir de una serie de variables.

Otro ejemplo podría ser predecir la renta de una persona a partir de features, que es el ejemplo que vamos a hacer. El caso es que en regresión la label (el valor a predecir) se le suele llamar target, y debe ser una variable continua (numérica, y a ser posible con decimales). Todos los conceptos de sobreajuste, generalización y similares aplican también a los problemas de regresión.

Sin más dilación, comencemos por generar un dataset sintético para regresión. sklearn contiene la función sklearn.datasets.make_regression, pero los datasets que genera son demasiado sencillos: suele ser muy fácil ajustar el problema de regresión. Así que vamos a usar álgebra con Numpy para generar un dataset más interesante para regresión.

En los problemas de clasificación, como hemos visto hasta ahora, se puede visualizar bien "3 dimensiones en 2D". Quiero decir: teníamos dos features y la label, lo cual en realidad suman 3 dimensiones distintas. Lo visualizábamos dibujando en 2D las features, y poniendo la label como color, a modo de la tercera dimensión.

En regresión suele ser más complicado visualizarlo así, así que vamos a generar solo una feature que irá en el eje x, el target que irá en el eje y, y el color de los puntos será uniforme e indiferente. Vamos a ello: pondremos en el eje x los años trabajando de distintos individuos, y en el eje y el target, que será el salario bruto anual:

In [69]:
# Comencemos por generar la coordenada x:
np.random.seed(12)
x = np.random.normal(loc=9, scale=16, size=3200)

# Vamos a quitar todos los valores que sean menores de 0
# (porque nadie puede llevar -1 años trabajados):
x = x[x>=0]

# Nos han quedado estas observaciones:
print("Número de observaciones: " + str(len(x)))

# Y ahora la coordenada y, es decir, el target,
# a modo de función de x. Y vamos a añadir 
# algo de ruido metiendo algo de random:
y = 8000*x**(1.0/3.0) + 0.8*x**3 + 21000 + np.random.normal(scale=5000, size=len(x))

# A DataFrame:
dataset_regresion = pd.DataFrame({"años trabajados": x, "salario bruto": y})

# Primeras cinco filas de nuestro dataset de regresión:
dataset_regresion[:5]
Número de observaciones: 2232
Out[69]:
años trabajados salario bruto
0 16.567773 49222.893302
1 12.879032 38385.332814
2 21.050285 46799.648090
3 9.082033 37556.618248
4 7.076357 36261.284893
In [70]:
dataset_regresion.plot.scatter(x="años trabajados",y="salario bruto",figsize=(8,6))
pass

Dividimos en train y test:

In [71]:
dataset_regresion_spliteado = train_test_split(dataset_regresion,
                                                   train_size=0.8, # 80% training
                                                   test_size=0.2,  # 20% testing
                                                   random_state=24 # seed para el generador de num. al.
                                                  )

dataset_regresion_training = dataset_regresion_spliteado[0]
dataset_regresion_testing = dataset_regresion_spliteado[1]

Y los dibujamos:

In [72]:
plt.figure(figsize=(16,7))

# Training
axes1=plt.subplot(1,2,1, title="Training")
dataset_regresion_training.plot.scatter(x="años trabajados",y="salario bruto",ax=axes1)

# Testing
axes2=plt.subplot(1,2,2, title="Testing")
dataset_regresion_testing.plot.scatter(x="años trabajados",y="salario bruto",ax=axes2)
pass

No es súper realista, pero nos servirá para nuestros propósitos.

Sklearn incluye muchos modelos regresores. De hecho, prácticamente todos los modelos clasificadores tienen su variante regresora. Estos son algunos de ellos:

Algoritmo Módulo Comentario
Mínimos Cuadrados linear_model.LinearRegression Regresión lineal con Least Squares
Ridge Regression linear_model.Ridge Least Squares con regularización $L2$
Lasso Regression linear_model.Lasso Least Squares con regularización $L1$
Elastic Net linear_model.ElasticNet Least Squares con regularización $L1$ y $L2$
KNN Regression neighbors.KNeighborsRegressor KNN adaptado a regresión (filtro de media)
Radius NN Regression neighbors.RadiusNeighborsRegressor NN teniendo en cuenta el radio y no el número de vecinos de la muestra
Regression Trees tree.DecisionTreeRegressor Árboles de decisión para regresión
Linear SVM Regresión svm.LinearSVR Implementación de SVM para regresión. Escala mejor para un gran número de muestras. Emplea la biblioteca liblinear
SVM Regression svm.SVR Implementación de SVM para regresión. Admite kernels. Emplea la biblioteca libsvm
Isotonic Regression isotonic.IsotonicRegression Regresión diseñada para datos con valor creciente
Procesos Gaussianos gaussian_process.GaussianProcess Método de regresión con salida probabilística
Regresion Polinomica preprocessing.PolynomialFeatures Expansión de los datos para regresión no lineal

Comencemos a hacer regresiones:

6.A. Regresión lineal simple (mínimos cuadrados ordinarios)

En el ejercicio de Numpy aprendimos a hacer con álgebra lineal el ajuste por mínimos cuadrados. Sklearn nos permite hacerlo de forma sencilla con sklearn.linear_model.LinearRegression:

In [73]:
from sklearn.linear_model import LinearRegression

ols = LinearRegression(fit_intercept=True, # Añadir intercept, bias u ordenada en el origen
                       normalize=False, # Normalizar las features antes de entrenar. Es el valor por
                       )                # defecto, pero lo pongo para que lo veas.

Exactamente igual que con los clasificadores, llamamos al método .fit() para entrenar la regresión:

In [74]:
ols.fit(X=dataset_regresion_training[["años trabajados"]], # Aunque sea una sola feature,
        y=dataset_regresion_training["salario bruto"])     # sklearn quiere que la seleccionemos
                                                           # como una lista de columnas, con solo
                                                           # un elemento
Out[74]:
LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)

El problema de regresión lineal planteado es el siguiente:

$$salario = \beta_{0} + \beta_{1}años$$

Ya hemos entrenado el modelo, así que tenemos los parámetros $\beta$ estimados:

In [75]:
print("beta0 o intercept: " + str(ols.intercept_))
print("beta1: " + str(ols.coef_[0]))
beta0 o intercept: 19842.4003843
beta1: 1799.7291852

El rendimiento en regresión no se puede medir en acertar/fallar el label, puesto que se trata de números continuos. Así, podemos medir cuánto nos hemos acercado con nuestra predicción al número real. Sklearn tiene en sus regresores (una vez entrenados) el método .score(), exactamente igual que en los clasificadores.

No obstante, cuando aplicamos .score(), sklearn calcula el $R^2$ (coeficiente de determinación), el cual NO recomiendo como medida del rendimiento de nuestra regresión. Vamos a probarlo de todos modos:

In [76]:
print("R cuadrado : " + str(ols.score(X=dataset_regresion_testing[["años trabajados"]],
                                      y=dataset_regresion_testing["salario bruto"])))
R cuadrado : 0.81113119225

En problemas de regresión, recomiendo utilizar otro tipo de métricas, como puede ser el error cuadrático medio (mean squared error o $MSE$), o bien la raíz de este valor (root mean squared error o $RMSE$). Afortunadamente, sklearn tiene varias métricas para medir el rendimiento de regresores en sklearn.metrics.

Vamos a calcular el $MSE$ y $RMSE$ utilizando sklearn.metrics.mean_squared_error:

In [77]:
# mean_squared_error toma como primer argumento los valores
# reales del target, y como segundo los valores que predecimos:
from sklearn.metrics import mean_squared_error

# Medimos el MSE en el conjunto de testing.
mse = mean_squared_error(y_true=dataset_regresion_testing["salario bruto"],
                         y_pred=ols.predict(dataset_regresion_testing[["años trabajados"]]))

# El RMSE es la raíz de este valor. Podemos utilizar
# np.sqrt o math.sqrt (de la biblioteca estándar)
# para calcularlo:
rmse = np.sqrt(mse)

print("MSE: " + str(mse))
print("RMSE: " +  str(rmse))
MSE: 113092578.256
RMSE: 10634.4994361

Y este es nuestro $E_{out}$ para nuestra regresión lineal. El $RMSE$ se puede entender como: en términos medios, fallamos las predicciones por unos 10.600€ de salario (aprox), ya sea de más o de menos, sobre el valor real.

También podemos dibujar nuestra línea de regresión:

In [78]:
figura = plt.figure(figsize=(16,7))
figura.suptitle("Regresión lineal simple (mínimos cuadrados ordinarios)", fontsize=14, fontweight="bold")

# Training
axes1=plt.subplot(1,2,1, title="Training")
dataset_regresion_training.plot.scatter(x="años trabajados",y="salario bruto",ax=axes1)
# Dibujamos la recta de regresión:
plt.plot(dataset_regresion_training["años trabajados"],
         ols.predict(dataset_regresion_training[["años trabajados"]]),
         c="red")

# Testing
axes2=plt.subplot(1,2,2, title="Testing")
dataset_regresion_testing.plot.scatter(x="años trabajados",y="salario bruto",ax=axes2)
# Dibujamos la recta de regresión:
plt.plot(dataset_regresion_testing["años trabajados"],
         ols.predict(dataset_regresion_testing[["años trabajados"]]),
         c="red")
pass

El ajuste no es ideal. La regresión lineal simple no puede hacer mucho más...

Teniendo en cuenta la "forma" de nuestro dataset, ningún modelo lineal va a poder darnos resultados mucho mejores. Dicho esto, vamos a mencionar brevemente tres de las regularizaciones más famosas para la regresión lineal:

6.B. Regularización: Ridge, LASSO y Elastic Net

En clasificación vimos que algunos modelos toman como argumento algunos parámetros como "l1" o "l2". Esto son lo que se llaman regularizaciones. Otros eran simplemente argumentos que especificaban las características del modelo (como puede ser el número de vecinos a contemplar en KNN, o la profundidad del árbol en los árboles de decisión).

La regularización es un proceso que intenta alterar ligeramente la formulación matemática de un modelo de machine learning, con la intención de prevenir el sobreajuste. Realmente, consiste en añadir hiperparámetros a nuestro modelo que son elegidos por nosotros, que incluyen ciertas penalizaciones. Dichas penalizaciones hacen que cuando el modelo se entrene, lo haga acorde a reducir el error (como siempre). Pero ahora este error ha sido alterado, incluyendo dichas penalizaciones; de forma que el modelo intentará conseguir un balance entre ajustar bien, y no ser penalizado en exceso.

Diseñar una regularización es un proceso algo complicado: requiere saber qué se quiere penalizar, cómo añadirlo a la formulación matemática del modelo, y medir los resultados. Dicho esto, existen tantas posibles regularizaciones como las que se nos pueda ocurrir.

Tanto para la regresión lineal como para la regresión logística (y también aplicable a otros varios modelos) existen una serie de regularizaciones que se conoce que funcionan bastante bien y están bastante consolidadas. Las más famosas y utilizadas son estas dos:

  1. Regularización o penalización $L2$: no entraré en terreno matemático. A efectos prácticos, la regularización $L2$ consigue que los parámetros estimados por nuestro modelo (las $\beta$ según la nomenclatura que hemos ido siguiendo) no tengan un valor absoluto demasiado grande. Es decir: introducir la penalización $L2$ en una regresión lineal hará que el modelo "sea forzado" a estimar parámetros $\beta$ pequeños (0.5, 3.1, -0.6) antes que parámetros $\beta$ grandes (3.320, 4.930, -14.291). Debido a esto, esta regularización se conoce también como parameter shrinking (contracción de parámetros estimados; puesto que los hace chiquititos) . El nombre "científico" es regularización de Tikhonov.
  2. Regularización o penalización $L1$: es más curiosa aún. Consigue que algunos de los parámetros $\beta$ estimados sean 0, es decir, que no sean útiles ni afecten a la predicción. Así que en el fondo, lo que hace es forzar la selección de features o variables dentro del modelo. Que yo sepa, no tiene nombre "científico".

Como estas dos regularizaciones son tan famosas, las regresiones lineales que las aplican tienen nombres específicos (aunque repito que se pueden utilizar en modelos más allá de las regresiones lineales):

  1. Regresión Ridge: es la regresión lineal que aplica la regularización $L2$.
  2. Regresión LASSO (o LASSO a secas): es la regresión lineal que aplica la regularización $L1$.
  3. Elastic Net: la regresión lineal que aplica tanto $L1$ como $L2$, y nos permite fijar qué proporción de cada una de ellas aplica.

A cada una de estas regularizaciones les acompaña por tanto un hiperparámetro que fijamos nosotros, y que dicta la "intensidad" de estas regularizaciones/penalizaciones.

Mucha teoría. Vamos a hacer un ejemplo utilizando la regresión Ridge, disponible en sklearn en sklearn.linear_model.Ridge:

In [79]:
from sklearn.linear_model import Ridge

# Formulamos la Ridge:
ridge = Ridge(alpha=1, # Hiperparámetro de la regularización, que marca la intensidad del shrinking
              fit_intercept=True)

# La entrenamos:
ridge.fit(X=dataset_regresion_training[["años trabajados"]],
        y=dataset_regresion_training["salario bruto"])
Out[79]:
Ridge(alpha=1, copy_X=True, fit_intercept=True, max_iter=None,
   normalize=False, random_state=None, solver='auto', tol=0.001)

Vamos a ver si los parámetros $\beta$ estimados efectivamente se han hecho más pequeños, comparados con los de la regresión lineal simple por mínimos cuadrados que hicimos antes. Prestamos atención a $\beta_1$ (y no al $\beta_0$ o intercept):

In [80]:
print("beta0 o intercept: " + str(ridge.intercept_))
print("beta1: " + str(ridge.coef_[0]))
beta0 o intercept: 19842.5314253
beta1: 1799.72122671

Algo lo ha reducido, pero no mucho. Si aumentamos el hiperparámetro alpha, es decir, la intensidad de la regularización:

In [81]:
# Formulamos la Ridge:
ridge2 = Ridge(alpha=100000, # No hagas esto nunca; es demsiado alto...
               fit_intercept=True)

# La entrenamos:
ridge2.fit(X=dataset_regresion_training[["años trabajados"]],
        y=dataset_regresion_training["salario bruto"])

print("beta0 o intercept: " + str(ridge2.intercept_))
print("beta1: " + str(ridge2.coef_[0]))
beta0 o intercept: 28928.5882408
beta1: 1247.89977467

Efectivamente, lo va reduciendo. En el caso de que tuviéramos más features (y por tanto más estimadores), el efecto sería mucho más notable. Cuando tenemos solo una variable, la regularización no surte tanto efecto; de forma que lo que hemos hecho aquí es poco útil, y solo tiene propósitos educativos.

Siéntete libre de probar también con LASSO y Elastic Net para ver cómo cambian los parámetros estimados, a ser posible con un problema multivariante (es decir, con más de una feature) para ver mejor los efectos de las regularizaciones.

La pregunta que queda es: ¿cómo elegir una regularización, así como la intensidad de la misma? La respuesta es la misma que cuando se trata de seleccionar modelos (al fin y al cabo, sigue siendo seleccionar el modelo, o su variante regularizada, que mejor funciona). Para ello, veremos luego el proceso de validación de modelos y su selección.

Ahora, continuamos con las regresiones. Vamos a probar la polinómica.

6.C. Regresión polinómica: transformación no lineal de nuestro dataset

Las regularizaciones sobre un modelo lineal siguen dándonos un modelo lineal. Vamos a probar a hacer regresiones no lineales utilizando polinomios. sklearn.preprocessing incluye las PolynomialFeatures, que básicamente coge nuestras features y las transforma en polinomios del grado que le pidamos.

Normalmente esta función es más útil cuando tenemos más de una feature. Supongamos que tenemos dos (años trabajados y altura de la persona):

$$salario = \beta_{0} + \beta_{1}años + \beta_{2}altura$$

Aplicando PolynomialFeatures de grado 2, nuestro problema se transformaría en: $$salario = \beta_{0} + \beta_{1}años + \beta_{2}altura + \beta_{3}años^2 + \beta_{4}altura^2 + \beta_{5}(años\times altura)$$

Lo cual nos podría dar un mejor ajuste, ya que contempla mayor nivel de interacción entre las features.

En nuestro caso donde solo tenemos una feature, la expansión polinómica de grado 2 generará:

$$salario = \beta_{0} + \beta_{1}años + \beta_{2}años^2$$

Vamos a probarlo:

In [82]:
from sklearn.preprocessing import PolynomialFeatures

# Ahora, creamos una instancia de PolynomialFeatures
# de grado 2:
poly2 = PolynomialFeatures(degree=2,
                           include_bias=False) # Ya añadiremos el bias/intercept cuando entrenemos modelo

# Y ahora transformamos nuestro dataset, utilizando
# un método para los transformadores de features
# llamado fit_transform():
features_poly2 = poly2.fit_transform(dataset_regresion_training[["años trabajados"]])

A partir de lo que solo era una feature (años trabajados), se han creado dos:

In [83]:
features_poly2
Out[83]:
array([[  18.32547077,  335.82287898],
       [   8.55005136,   73.10337821],
       [  14.04873749,  197.36702519],
       ..., 
       [   6.10715458,   37.29733706],
       [  29.69666446,  881.89188008],
       [   1.22215377,    1.49365984]])

Vamos a meterlas en un DataFrame nuevo:

In [84]:
dataset_poly2_training = pd.DataFrame(features_poly2)

# Nombramos columnas:
dataset_poly2_training.columns = ["años trabajados", "años trabajados ^2"]

# Y añadimos los targets del conjunto training original, inalterados.
# OJO: queremos que pandas "pase" de los índices porque nos desordenaría
# la columna de los targets con respecto al resto, así que
# quitamos los índices pasando la columna a array de Numpy:
dataset_poly2_training["salario bruto"] = np.array(dataset_regresion_training["salario bruto"])

Así era el conjunto de training antes de las PolynomialFeatures:

In [85]:
dataset_regresion_training[:7]
Out[85]:
años trabajados salario bruto
1691 18.325471 49640.369885
1020 8.550051 43213.503246
734 14.048737 44235.170764
1816 16.467413 48466.798982
517 16.806036 44766.191453
1172 27.555046 66515.452648
706 2.922094 34130.001555

Y así es después:

In [86]:
dataset_poly2_training[:7]
Out[86]:
años trabajados años trabajados ^2 salario bruto
0 18.325471 335.822879 49640.369885
1 8.550051 73.103378 43213.503246
2 14.048737 197.367025 44235.170764
3 16.467413 271.175683 48466.798982
4 16.806036 282.442856 44766.191453
5 27.555046 759.280568 66515.452648
6 2.922094 8.538632 34130.001555

Muy bien. Ahora entrenemos la regresión. El trabajo de convertir la regresión lineal en regresión polinómica ya está hecho: ya hemos aplicado la transformación no lineal a nuestra feature:

Aplicar una regresión (o modelo) lineal en un espacio dimensional (o dataset) resultado de la transformación no lineal de las features, da como resultado un modelo no lineal.

Entrenamos (regresión lineal simple de mínimos cuadrados ordinarios):

In [87]:
lr_2 = LinearRegression(fit_intercept=True)

lr_2.fit(X=dataset_poly2_training[["años trabajados", "años trabajados ^2"]],
         y=dataset_poly2_training["salario bruto"])
Out[87]:
LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)

Hecho. Ahora, no obstante, tenemos también que aplicar la transformación al conjunto de testing, ya que para predecir, tenemos igualmente que realizar la transformación polinómica a las features:

In [88]:
# Hacemos fit_transform del PolynomialFeatures
# y metemos en DataFrame nuevo:
dataset_poly2_testing = pd.DataFrame(poly2.fit_transform(dataset_regresion_testing[["años trabajados"]]))

# Renombramos columnas:
dataset_poly2_testing.columns = ["años trabajados", "años trabajados ^2"]

# Y añadimos los targets:
dataset_poly2_testing["salario bruto"] = np.array(dataset_regresion_testing["salario bruto"])

Conjunto de testing antes de la transformación:

In [89]:
dataset_regresion_testing[:7]
Out[89]:
años trabajados salario bruto
924 17.448235 43638.733378
2116 12.033330 39952.371043
609 29.952618 70986.589809
906 7.590485 38425.022364
529 9.577510 35740.035714
664 3.292486 24966.195502
92 31.527176 66513.641263

Y después:

In [90]:
dataset_poly2_testing[:7]
Out[90]:
años trabajados años trabajados ^2 salario bruto
0 17.448235 304.440914 43638.733378
1 12.033330 144.801027 39952.371043
2 29.952618 897.159344 70986.589809
3 7.590485 57.615469 38425.022364
4 9.577510 91.728705 35740.035714
5 3.292486 10.840464 24966.195502
6 31.527176 993.962810 66513.641263

Medimos el $RMSE$ en el conjunto de testing:

In [91]:
rmse_poly2 = np.sqrt(mean_squared_error(y_true=dataset_poly2_testing["salario bruto"],
                                        y_pred=lr_2.predict(dataset_poly2_testing[["años trabajados","años trabajados ^2"]])))
print("RMSE: " +  str(rmse_poly2))
RMSE: 5660.35464757

Parece que hemos conseguido dejar el $E_{out}$ a la mitad con nuestra regresión polinómica (comparando con la regresión lineal simple). Dibujemos:

In [92]:
# Cuidado: para dibujar curvas, es importante
# que los datos en el eje x estén ordenados de menor a mayor
# (y por tanto, re-ordenar también los del)
# eje y acorde. En este caso, como las coordenadas
# de nuestras curvas se obtienen con lr_2.predict(),
# con ordenar los datos del eje x (años trabajados) basta:

figura = plt.figure(figsize=(16,7))
figura.suptitle("Regresión polinómica (transformación polinómica de las features)", fontsize=14, fontweight="bold")


### Training
axes1=plt.subplot(1,2,1, title="Training")

# Dibujamos los puntos (aquí el orden da igual):
dataset_poly2_training.plot.scatter(x="años trabajados",y="salario bruto",ax=axes1)

# Ordenamos los datos por el eje x con el método
# .sort_values() de los DataFrames:
dataset_poly2_training_sorted = dataset_poly2_training.sort_values(by="años trabajados")

# Dibujamos la curva de regresión polinómica:
plt.plot(dataset_poly2_training_sorted["años trabajados"],
         lr_2.predict(dataset_poly2_training_sorted[["años trabajados","años trabajados ^2"]]),
         c="red")


### Testing
axes2=plt.subplot(1,2,2, title="Testing")

# Dibujamos los puntos
dataset_poly2_testing.plot.scatter(x="años trabajados",y="salario bruto",ax=axes2)

# Y también ordenamos los de testing:
dataset_poly2_testing_sorted = dataset_poly2_testing.sort_values(by="años trabajados")

# Dibujamos la curva de regresión polinómica:
plt.plot(dataset_poly2_testing_sorted["años trabajados"],
         lr_2.predict(dataset_poly2_testing_sorted[["años trabajados","años trabajados ^2"]]),
         c="red")
pass

No dibujamos el cuadrado de los años trabajados. Para temas de visualización, queda mejor si solo utilizamos las features originales (realmente, hemos usado un modelo lineal sobre un dataset transformado de forma no lineal).

Vemos que ajusta mucho mejor. Podríamos probar también con polinomios de grado 3, 4, 5... Pero si elevamos mucho el grado del polinomio, acabaríamos sobreajustando el modelo de nuevo.

Probemos otros regresores distintos.

6.D. Otros regresores

Sklearn incluye más regresores aparte de la regresión lineal, su posible expansión polinómica y sus regularizaciones. Veamos algunos:

6.D.a. Árboles regresores

sklearn.tree.DecisionTreeRegressor nos permite usar árboles para problemas de regresión. Los árboles son modelos no lineales, así que podemos olvidarnos de tema de transformaciones no lineales del dataset:

In [93]:
from sklearn.tree import DecisionTreeRegressor

# Formulamos árbol:
tree_reg = DecisionTreeRegressor(criterion="mse", # Al ser árboles regresores, no usan gini o entropía, sino MSE
                                 max_depth=4,
                                 min_samples_split=5,
                                 min_samples_leaf=5
                                )

# Entrenamos:
tree_reg.fit(X=dataset_regresion_training[["años trabajados"]],
             y=dataset_regresion_training["salario bruto"])

# Medimos:
mse = mean_squared_error(y_true=dataset_regresion_testing["salario bruto"],
                         y_pred=tree_reg.predict(dataset_regresion_testing[["años trabajados"]]))
rmse = np.sqrt(mse)
print("MSE: " + str(mse))
print("RMSE: " +  str(rmse))
MSE: 31610245.4297
RMSE: 5622.29894525

Dibujamos:

In [94]:
figura = plt.figure(figsize=(16,7))
figura.suptitle("Árbol regresor (profundidad máxima de 4)", fontsize=14, fontweight="bold")

# De nuevo, al no ser una línea recta
# la solución, debemos utilizar el dataset
# ordenado en el eje x de menor a mayor
# (ya lo tenemos ordenado de antes):

### Training
axes1=plt.subplot(1,2,1, title="Training")

# Dibujamos los puntos
dataset_poly2_training.plot.scatter(x="años trabajados",y="salario bruto",ax=axes1)

# Dibujamos la regresión de árbol (con el dataset ordenador por la x):
plt.plot(dataset_poly2_training_sorted["años trabajados"],
         tree_reg.predict(dataset_poly2_training_sorted[["años trabajados"]]),
         c="red")

### Testing
axes2=plt.subplot(1,2,2, title="Testing")

# Dibujamos los puntos
dataset_poly2_testing.plot.scatter(x="años trabajados",y="salario bruto",ax=axes2)

# Dibujamos la regresión de árbol:
plt.plot(dataset_poly2_testing_sorted["años trabajados"],
         tree_reg.predict(dataset_poly2_testing_sorted[["años trabajados"]]),
         c="red")
pass

Las reglas de decisión y splits de este árbol se pueden dibujar de la misma forma que las de los árboles clasificadores (cualquier modelo que proviene de sklearn.tree puede ser dibujado con el mismo procedimiento).

6.D.b. Ensembles

Al igual que existen ensembles de clasificadores, existen ensembles de regresores. Para los ensembles de bagging la técnica suele ser hacer la media de las predicciones de cada modelo lanzado. En boosting, se intenta ir corrigiendo el $MSE$ de los modelos anteriores con los nuevos.

Probemos un Random Forest regresor sklearn.ensemble.RandomForestRegressor:

In [95]:
from sklearn.ensemble import RandomForestRegressor

rf_reg = RandomForestRegressor(n_estimators=20, # 20 árboles
                               criterion="mse",
                               max_depth=3,
                               min_samples_split=10,
                               min_samples_leaf=5,
                               bootstrap=True)

# Entrenamos:
rf_reg.fit(X=dataset_regresion_training[["años trabajados"]],
             y=dataset_regresion_training["salario bruto"])

# Medimos:
mse = mean_squared_error(y_true=dataset_regresion_testing["salario bruto"],
                         y_pred=rf_reg.predict(dataset_regresion_testing[["años trabajados"]]))
rmse = np.sqrt(mse)
print("MSE: " + str(mse))
print("RMSE: " +  str(rmse))
MSE: 28110566.4526
RMSE: 5301.93987636

Dibujamos:

In [96]:
figura = plt.figure(figsize=(16,7))
figura.suptitle("Random Forests (regresor)", fontsize=14, fontweight="bold")

### Training
axes1=plt.subplot(1,2,1, title="Training")

# Dibujamos los puntos
dataset_poly2_training.plot.scatter(x="años trabajados",y="salario bruto",ax=axes1)

# Dibujamos la regresión:
plt.plot(dataset_poly2_training_sorted["años trabajados"],
         rf_reg.predict(dataset_poly2_training_sorted[["años trabajados"]]),
         c="red")

### Testing
axes2=plt.subplot(1,2,2, title="Testing")

# Dibujamos los puntos
dataset_poly2_testing.plot.scatter(x="años trabajados",y="salario bruto",ax=axes2)

# Dibujamos la regresión:
plt.plot(dataset_poly2_testing_sorted["años trabajados"],
         rf_reg.predict(dataset_poly2_testing_sorted[["años trabajados"]]),
         c="red")
pass

6.D.c. Regresión por Nearest Neighbors

Su funcionamiento es sencillo: mira los vecinos cercanos a los puntos a predecir, y hace la media de los valores de dichos vecinos. Al igual que su homólogo clasificador, es un modelo no lineal.

Se encuentra en sklearn.neighbors.KNeighborsRegressor:

In [97]:
from sklearn.neighbors import KNeighborsRegressor

knn_reg = KNeighborsRegressor(n_neighbors=5, # 5-NN
                              metric="euclidean") 

# Entrenamos:
knn_reg.fit(X=dataset_regresion_training[["años trabajados"]],
             y=dataset_regresion_training["salario bruto"])

# Medimos:
mse = mean_squared_error(y_true=dataset_regresion_testing["salario bruto"],
                         y_pred=knn_reg.predict(dataset_regresion_testing[["años trabajados"]]))
            
rmse = np.sqrt(mse)
print("MSE: " + str(mse))
print("RMSE: " +  str(rmse))
MSE: 25946008.8913
RMSE: 5093.72249846
In [98]:
figura = plt.figure(figsize=(16,7))
figura.suptitle("5-NN (regresor)", fontsize=14, fontweight="bold")

### Training
axes1=plt.subplot(1,2,1, title="Training")

# Dibujamos los puntos
dataset_poly2_training.plot.scatter(x="años trabajados",y="salario bruto",ax=axes1)

# Dibujamos la regresión:
plt.plot(dataset_poly2_training_sorted["años trabajados"],
         knn_reg.predict(dataset_poly2_training_sorted[["años trabajados"]]),
         c="red") # En rojo se ve mejor

### Testing
axes2=plt.subplot(1,2,2, title="Testing")

# Dibujamos los puntos
dataset_poly2_testing.plot.scatter(x="años trabajados",y="salario bruto",ax=axes2)

# Dibujamos la regresión:
plt.plot(dataset_poly2_testing_sorted["años trabajados"],
         knn_reg.predict(dataset_poly2_testing_sorted[["años trabajados"]]),
         c="red") # En rojo se ve mejor
pass

6.D.d. Support Vector Regression

Aunque no tan famosos como sus homólogos clasificadores, los SVRs (Support Vector Regressors) también funcionan y generalizan muy bien. El que vamos a usar está en sklearn.svm.SVR:

In [99]:
from sklearn.svm import SVR

svr = SVR(C=0.5,
          kernel="poly", # kernel polinómico
          degree=2,      # de segundo grado
         )

# Entrenamos:
svr.fit(X=dataset_regresion_training[["años trabajados"]],
             y=dataset_regresion_training["salario bruto"])

# Medimos:
mse = mean_squared_error(y_true=dataset_regresion_testing["salario bruto"],
                         y_pred=svr.predict(dataset_regresion_testing[["años trabajados"]]))
            
rmse = np.sqrt(mse)
print("MSE: " + str(mse))
print("RMSE: " +  str(rmse))
MSE: 34760292.7282
RMSE: 5895.78601445
In [100]:
figura = plt.figure(figsize=(16,7))
figura.suptitle("SVR (kernel polinómico de grado 2)", fontsize=14, fontweight="bold")

### Training
axes1=plt.subplot(1,2,1, title="Training")

# Dibujamos los puntos
dataset_poly2_training.plot.scatter(x="años trabajados",y="salario bruto",ax=axes1)

# Dibujamos la regresión:
plt.plot(dataset_poly2_training_sorted["años trabajados"],
         svr.predict(dataset_poly2_training_sorted[["años trabajados"]]),
         c="red") # En rojo se ve mejor

### Testing
axes2=plt.subplot(1,2,2, title="Testing")

# Dibujamos los puntos
dataset_poly2_testing.plot.scatter(x="años trabajados",y="salario bruto",ax=axes2)

# Dibujamos la regresión:
plt.plot(dataset_poly2_testing_sorted["años trabajados"],
         svr.predict(dataset_poly2_testing_sorted[["años trabajados"]]),
         c="red") # En rojo se ve mejor
pass

Hay muchos más. Cada vez que veas un modelo con la coletilla Regressor... Ya sabes :)

7. Clustering

A diferencia de la clasificación y la regresión, los modelos de clustering entran dentro de lo que se llama aprendizaje no supervisado. Tanto en clasificación como en regresión, entrenamos los modelos con datos de los que sabemos el valor real a predecir:

  • En el caso de clasificación los hemos llamado labels o clases (calvo/no-calvo en nuestro ejemplo).
  • En el caso de regresión los hemos llamado targets (el salario bruto anual).

En el aprendizaje no supervisado no es así: no tenemos un valor real a predecir: bien porque no está disponible en los datos, o bien porque no existe como tal. Por eso en aprendizaje no supervisado no se suele hablar de modelos precisos o imprecisos, sino de modelos útiles o apropiados para lo que queremos conseguir.

Un ejemplo de aprendizaje no supervisado son los modelos de clustering. Su principio es sencillo: buscar entre nuestros datos, y agruparlos de forma que cada grupo sea lo más homogéneo posible, y lo más diferente posible del resto de grupos.

A mí me gusta dividir los algoritmos/modelos de clustering en dos grandes grupos:

  1. Clusterings basados en distancia: al fin y al cabo, nuestras observaciones se pueden representar como puntos en un espacio dimensional, donde las dimensiones son tantas como el número de features que tenemos. Los modelos de clustering basados en distancia utilizan una métrica de distancia (como puede ser la euclídea, la Manhattan o la Mahalanobis) para separar los distintos grupos o clusters.
  2. Clusterings basados en densidad: computan los grupos no teniendo en cuenta solo la distancia absoluta entre los puntos, sino detectando "zonas" donde la densidad de puntos/observaciones es mayor.

Sklearn tiene modelos de clustering de ambos tipos. Preparemos un par de datasets, y pongámonos manos a la obra:

7.A. Los datos

Utilizaremos dos datasets: por un lado uno sintético creado con sklearn.datasets.make_blobs(); y por otro lado otro preparado para probar algoritmos de clustering. Comencemos por el de make_blobs:

In [101]:
dataset = make_blobs(n_samples=1350,
                     n_features=2,
                     centers=6,
                     cluster_std=1.5,
                     random_state=11)

Recuerdo que make_blobs nos devuelve una tupla con dos arrays de Numpy: uno con las features, y otro con las labels. Vamos a ignorar completamente las labels, al estar en aprendizaje no supervisado. Vamos a cargar simplemente las features en un DataFrame:

In [102]:
dataset_clust = pd.DataFrame(dataset[0])

# Nombramos columnas. Paso de inventarme historias:
# les llamo x e y a las features:
dataset_clust.columns = ["x","y"]

dataset_clust[:5]
Out[102]:
x y
0 -8.452267 -9.652973
1 -2.003873 5.111289
2 8.914786 6.562965
3 -11.982325 0.555952
4 -7.190894 -10.915288

Dibujamos:

In [103]:
dataset_clust.plot(kind="scatter",
                     x="x",
                     y="y",
                     figsize=(8,6))
pass

El random_state o seed del generador de números aleatorios ha sido elegido a propósito, para que los dos blobs del centro parezcan más o menos uno solo. Pero no me gusta: demasiado sencillo.

Vamos a generar algo de ruido, metiendo un blob grande y con gran desviación típica, y lo añadimos al dataset, a modo de "datos huérfanos", sin cluster ovbio:

In [104]:
ruido = make_blobs(n_samples=550,
                   n_features=2,
                   centers=1,
                   cluster_std=8,
                   center_box=[0,0],
                   random_state=19)

# Lo metemos en un DataFrame:
ruido_df = pd.DataFrame(ruido[0])
ruido_df.columns = ["x","y"]

# Y lo concatenamos con dataset_clust:
dataset_clust_1 = pd.concat([dataset_clust, ruido_df])

Y dibujamos:

In [105]:
dataset_clust_1.plot(kind="scatter",
                     x="x",
                     y="y",
                     figsize=(8,6))
pass

Mucho mejor. Precioso.

Por otro lado, cargamos el dataset preparado (que resulta estar en formato csv). Estos datos son también sintéticos, cedidos amablemente por los creadores de la implementación en Python del HDBSCAN:

In [106]:
dataset_clust_2 = pd.read_csv("clustering_data.csv", sep=";")

dataset_clust_2[:5]
Out[106]:
x y
0 -0.121535 -0.228763
1 -0.220937 -0.252511
2 0.125904 -0.273143
3 -0.164537 -0.222244
4 -0.180824 -0.211075

Y lo dibujamos:

In [107]:
dataset_clust_2.plot(kind="scatter",
                     x="x",
                     y="y",
                     figsize=(8,6))
pass

Mola porque tiene clusters no esféricos, así como ruido.

Ya podemos empezar a lanzar modelos.

7.B. Algoritmos de clustering

Aquí muestro algunos de los que incluye sklearn (como digo siempre, hay más de los listados aquí):

Algoritmo Módulo
K-Means cluster.KMeans
Affinity propagation cluster.AffinityPropagation
Mean-shift cluster.MeanShif
Spectral clustering cluster.SpectralClustering
Agglomerative clustering (jerárquico) cluster.AgglomerativeClustering
DBSCAN cluster.DBSCAN
Gaussian mixtures mixture.GMM
Birch cluster.Birch

Probemos algunos de estos con nuestros dos datasets.

7.B.a. K-Means

El K-Means es el algoritmo de clustering más conocido. Sencillo de entender, basado en distancia euclídea y computacionalmente no muy intensivo (si bien se trata de un problema que se resuelve a través de métodos heurísticos, el Algoritmo de Lloyd hace un buen trabajo convergiendo para el K-Means en la mayoría de los casos). A lo largo de estas iteraciones, los centroides de los clusters se van moviendo, siendo la media de los puntos de cada cluster en la iteración anterior. Y de ahí el nombre K-means (o K-medias).

Usémoslo en nuestros dos datasets con sklearn.cluster.KMeans:

In [108]:
from sklearn.cluster import KMeans

km1 = KMeans(n_clusters=6, # Pega del KMeans: nos obliga a elegir número de clusters.
            max_iter=400)

# Entrenamos con el dataset_clust_1:
km1.fit(dataset_clust_1) # OJO: cada vez que se entrena, los clusters "no coinciden",
                         # es decir: lo que en un entrenamiento es el cluster número 3,
                         # en otro entrenamiento el cluster será el mismo (que contiene
                         # más o menos los mismos puntos), pero igual pasa a ser el número
                         # 2, o 4, o el que sea...
                         # Ocurre en todos los modelos de clustering no-determinísticos.

# Y predecimos a qué cluster pertenece cada punto
# que se representan con números, comenzando por el 0
# hasta n_clusters - 1:
result_kmeans_dataset1 = km1.predict(dataset_clust_1)
result_kmeans_dataset1
Out[108]:
array([4, 5, 0, ..., 4, 3, 3], dtype=int32)

Podemos utilizar estas etiquetas para pintar a qué cluster pertenece cada punto:

In [109]:
dataset_clust_1.plot(kind="scatter",
                     x="x",
                     y="y",
                     c=km1.predict(dataset_clust_1),
                     cmap=color_maps.Accent,
                     figsize=(9,7))
pass

Y para el otro dataset (donde asumimos que también hay 6 clusters):

In [110]:
km2 = KMeans(n_clusters=6,
             max_iter=600)

# Entrenamos con el dataset_clust_2:
km2.fit(dataset_clust_2)

# Y dibujamos:
dataset_clust_2.plot(kind="scatter",
                     x="x",
                     y="y",
                     c=km2.predict(dataset_clust_2),
                     cmap=color_maps.Accent,
                     figsize=(9,7))
pass

Rápidamente podemos ver los problemas del K-Means:

  1. Necesitamos decirle el número de clusters a crear. Existen técnicas para intentar inferir el número óptimo de clusters, pero en la práctica ninguna es maravillosa. Cierto es que en ocasiones el número de clusters es una imposición (por ejemplo, en una campaña de segmentación de marketing donde hay presupuesto para tratar de forma diferenciada a X grupos de clientes), pero en el resto de casos siempre es complicado...
  2. Los clusters generados siempre son esféricos (o circulares en el caso de dos dimensiones).
  3. Todos los puntos pertenecerán a algún cluster: el modelo no descarta ningún punto; por muy alejado que esté del resto. En casos muy extremos, puede hacer que un punto (o unos cuantos) muy alejados formen su propio cluster alejado; pero ningún punto es descartado.

7.B.b. Gaussian Mixture Model

El Gaussian Mixture Model (abreviado como GMM) es un modelo probabilístico que es usado comúnmente como modelo de clustering. De hecho, se trata de la generalización del K-means, o dicho de otro modo: el K-Means es un caso muy específico de GMM.

Probemos el sklearn.mixture.GaussianMixture (sí: no está en sklearn.cluster, sino en sklearn.mixture):

In [111]:
from sklearn.mixture import GaussianMixture

gmm1 = GaussianMixture(n_components=6, # En GMMs se habla de componentes en vez de clusters
                       covariance_type="diag", # Tipo de parámetros de covarianza a usar
                       max_iter=250) # Número máximo de iteraciones

# Entrenamos con el dataset_clust_1:
gmm1.fit(dataset_clust_1)

# Y dibujamos:
dataset_clust_1.plot(kind="scatter",
                     x="x",
                     y="y",
                     c=gmm1.predict(dataset_clust_1),
                     cmap=color_maps.Accent_r,
                     figsize=(9,7))
pass

GMM ha determinado el ruido como un cluster... A ver con el otro dataset:

In [112]:
gmm2 = GaussianMixture(n_components=6,
                       covariance_type="diag",
                       max_iter=500)

# Entrenamos con el dataset_clust_2:
gmm2.fit(dataset_clust_2)

# Y dibujamos:
dataset_clust_2.plot(kind="scatter",
                     x="x",
                     y="y",
                     c=gmm2.predict(dataset_clust_2),
                     cmap=color_maps.Accent,
                     figsize=(9,7))
pass

Bastante similar al K-Means en este caso. Consideraciones del GMM:

  1. Necesitamos definir el número de clusters a crear, igual que en el K-Means.
  2. Los clustes no tienen por qué ser esféricos: el K-Means contempla únicamente la media de los puntos para determinar los clusters (el vector de medias), mientras que el GMM contempla también la matriz de varianzas y covarianzas de los datos. Es esto lo que le permite la creación de clusters no esféricos.
  3. Todos los puntos pertenencen a algún cluster; no se descartan puntos.
  4. A diferencia del K-Means, la salida real del modelo no es una etiqueta de a qué componente/cluster pertenece cada punto; sino una distribución de probabilidad, donde se muestra la probabilidad de que cada punto pertenezca a cada componente/cluster. Eso sí: automáticamente, sklearn (y todas las bibliotecas que implementan GMM) nos da el componente/cluster al que es más probable que pertenezca el punto. No obstante, el método .predict_proba() nos ofrece esta distribución de probabilidad (¡otros modelos, incluso clasificadores, también lo ofrecen!)

7.B.c. Affinity Propagation

El Affinity Propagation es un modelo que se basa en grafos y que hace a cada punto votar a su cluster o exemplar favorito. A diferencia de en el K-Means y GMM, los centroides (o centros de los clusters) en el Affinity Propagation son puntos reales del propio dataset, los cuales son llamados exemplars. A partir de ahí, funciona de forma bastante similar al K-Means. Eso sí: la formulación del algoritmo permite utilizarlo con métricas no distanciales, esto es, que la métrica de similaridad/no-similaridad utilizada no tiene por qué satisfacer las cualidades de las distancias.

Utilicemos pues sklearn.cluster.AffinityPropagation:

In [113]:
from sklearn.cluster import AffinityPropagation

ap1 = AffinityPropagation(damping=0.9,  # Damping factor, >=0.5 y <1. A más grande, menos exemplars.
                          affinity="euclidean", # Métrica de distancia euclídea
                          preference=-1200.0, # Preferencias para cada punto (si es array), 
                                              # o generales(si es float). Por lo general,
                                              # recomiendo un número negativo, más grande o
                                              # más pequeño dependiendo de los datos. Si no,
                                              # nos sacará muchísimos clusters...
                          max_iter=3000) # Máximo de iteraciones

# Entrenamos con el dataset_clust_1:
ap1.fit(dataset_clust_1)

# Y dibujamos:
dataset_clust_1.plot(kind="scatter",
                     x="x",
                     y="y",
                     c=ap1.predict(dataset_clust_1),
                     cmap=color_maps.Accent,
                     figsize=(9,7))
pass

Sobre el otro dataset:

In [114]:
ap2 = AffinityPropagation(damping=0.95, 
                          affinity="euclidean",
                          preference=-3.0,
                          max_iter=3000)

# Entrenamos con el dataset_clust_2:
ap2.fit(dataset_clust_2)

# Y dibujamos:
dataset_clust_2.plot(kind="scatter",
                     x="x",
                     y="y",
                     c=ap2.predict(dataset_clust_2),
                     cmap=color_maps.Accent,
                     figsize=(9,7))
pass

Sobre el Affinity Propagation:

  1. No necesitamos decirle el número de clusters o exemplars: el propio modelo los determina a partir de la medida de preferencia y el dampening factor.
  2. Asume clusters esféricos.
  3. Todos los puntos pertenecen a algún cluster; no se descartan puntos.
  4. Es lento: tiene que determinar qué puntos van a ser exemplars, tanto en cantidad como en localización, y luego ya determinar a qué exemplar pertenece cada uno del resto de puntos.

7.B.d. Mean Shift

El Mean Shift es un modelo de clustering bastante potente. Básicamente el modelo asume que existe cierta función de densidad (probabilística) que explica la distribución de los datos, e intenta colocar centroides para maximizar la verosimilitud en función de los datos. Realiza dicha aproximación utilizando técnicas de KDE (Kernel Density Estimation).

El modelo se encuentra en sklearn.cluster.MeanShift:

In [115]:
from sklearn.cluster import MeanShift

ms1 = MeanShift(bandwidth=5.5, # El "ancho de banda" del kernel utilizado para la estimación
                cluster_all=False) # Podemos elegir que "clusterice" todos los puntos, o no

## Para los modelos de clustering sobre todo:
# existe un método llamado f.it_predict(), que
# entrena y predice sobre los datos dados automáticamente.
# No es ideal cuando tenemos training y testing, pero sí
# en modelos no supervisados. En vez de hacer .fit() y 
# luego .predict() sobre los mismos datos, hacemos
# .fit_predict() en un solo método:
dataset_clust_1.plot(kind="scatter",
                     x="x",
                     y="y",
                     c=ms1.fit_predict(dataset_clust_1),
                     cmap=color_maps.Accent_r, # Con el _r revertimos el orden de colores,
                     figsize=(9,7))            # lo cual en este caso nos pinta los -1 de color gris oscuro.
pass

¿Ves los puntos de color gris oscuro? Son los puntos descartados, que no pertenecen a ningún cluster (son etiquietados con el número $-1$ a modo de label en el resultado). Veámoslo en el otro dataset:

In [116]:
ms2 = MeanShift(bandwidth=0.17, 
                cluster_all=False)

# Entrenamos y dibujamos:
dataset_clust_2.plot(kind="scatter",
                     x="x",
                     y="y",
                     c=ms2.fit_predict(X=dataset_clust_2),
                     cmap=color_maps.Accent_r,
                     figsize=(9,7))
pass

Misma historia: los puntos gris oscuro han sido descartados.

Con el Mean Shift:

  1. No tenemos que determinar el número de clusters. Pero al igual que con el Affinity Propagation, tenemos que "tener cuidado" con los parámetros (el bandwidth en este caso, principalmente) para obtener una solución razonable.
  2. Los clusters resultantes son esféricos.
  3. Es capaz de descartar puntos, que considera que no pertenecen a ningún cluster (de forma opcional).

7.B.e. DBSCAN

El DBSCAN es un algoritmo de clustering basado en densidad: asume que hay clusters en regiones donde la densidad de los puntos es elevada. DBSCAN genera una transformación del espacio dimensional basándose en la densidad de los datos: los puntos en zonas densas son separados, mientras que los de zonas menos densas son apartados de forma iterativa.

Utilicemos sklearn.cluster.DBSCAN:

In [117]:
from sklearn.cluster import DBSCAN

dbs1 = DBSCAN(eps=0.95, # Parámetro epsilon, que cambia la "consideración de densidad" del algoritmo
              min_samples=15) # Número de puntos vecinos necesarios para que un punto sea "core" de densidad

# Entrenamos y dibujamos:
dataset_clust_1.plot(kind="scatter",
                     x="x",
                     y="y",
                     c=dbs1.fit_predict(dataset_clust_1),
                     cmap=color_maps.Accent_r, 
                     figsize=(9,7))         
pass

De nuevo: puntos grises, descartados. Probemos con el dataset complicado:

In [118]:
dbs2 = DBSCAN(eps=0.035,
              min_samples=20) 

# Entrenamos y dibujamos:
dataset_clust_2.plot(kind="scatter",
                     x="x",
                     y="y",
                     c=dbs2.fit_predict(dataset_clust_2),
                     cmap=color_maps.Accent_r, 
                     figsize=(9,7))        
pass

Sobre el DBSCAN:

  1. "Tocando" los parámetros epsilon y el número mínimo de puntos por core de densidad (no deberíamos hablar de centroides en DBSCAN), el modelo inferirá el número de clusters necesarios.
  2. Los clusters resultantes no tienen que ser esféricos (a la vista está).
  3. Descarta zonas de puntos donde la densidad es baja (descarta puntos por tanto).

Recuerdo: en el mundo del aprendizaje no supervisado (dentro del cual se incluye el clustering) NO hay una solución correcta; sino soluciones más convenientes que otras para un determinado propósito.

Recomiendo también ver (aparte del resto de modelos de clustering de sklearn) HDBSCAN, que es un modelo muy similar al DBSCAN. La diferencia es que HDBSCAN permite tener clusters densos, pero de densidades diferentes entre ellos.

No está incluído en sklearn, sino como biblioteca aparte. Sin embargo, se utiliza exactamente igual que los modelos de clustering de sklearn. E instalarlo es tan sencillo como:

$ pip install hdbscan

Nota a futuro para ti para siempre: sé que es tentador. Sé que te lo van a pedir en el trabajo. Sé que es... incluso divertido. Pero por favor: NO INCLUYAS NUNCA VARIABLES CATEGÓRICAS/DICOTÓMICAS EN UN MODELO DE CLUSTERING. La razón es sencilla; pongámonos en el caso de variables dicotómicas (que son o $1$ o $0$): en un dataset compuesto por dos features dicotómicas, de forma que cada observación puede ser $0|0$, $0|1$, $1|0$ o $1|1$... ¿Cuántos clusters hay? Pues 4, evidentemente: las 4 combinaciones. Introducir variables no numéricas en un modelo de clustering solo consigue una cosa: contaminar el modelo. Si el objetivo es hacer clustering y tienes variables numéricas y categóricas/dicotómicas, haz primero clustering sobre las numéricas, y luego computa todas las posibles combinaciones de variables categóricas/dicotómicas para cada cluster generado. Esto funciona, y es lo razonable.

Y con esto concluimos la parte de clustering. Seguimos con el feature engineering.

8. Feature Engineering

Lanzar modelos y obtener resultados es la parte sencilla en machine learning. El verdadero reto muchas veces es obtener features significativas y útiles para nuestros modelos.

Cuando te enfrentas a un problema de data science de primeras, normalmente comienzas por obtener los datos "en bruto": datos transaccionales mes a mes, comportamientos de cliente, imágenes de distintos tamaños y con diferentes perspectivas, etcétera.

Muchas veces el primer paso es generar variables derivadas a partir de estos datos en bruto (diferencias mes a mes, combinación de variables de crédito/débito, valores absolutos de balances, normalización del perfil RGB de fotografías...). Este tipo de transformaciones de los datos están muy vinculadas al sector/negocio en el que nos movemos, más que al punto de vista matemático.

Luego, es frecuente "filtrar" features: aquellas que tienen gran porcentaje de nulos, elevada correlación con otras, o valores constantes (que no permiten discriminar de ninguna manera).

A partir de ahí, hay tres formas principales de preparar los datasets para nuestros modelos:

  1. Incluir todas las variables (tanto en crudo o derivadas) válidas, y utilizar alguno de los modelos que "solo prestan atención" a las variables significativas. Ejemplos de estos modelos son aquellos que utilizan regularización $L1$ (como la regresión LASSO), los árboles (o ensembles de los mismos) y las redes.
  2. Utilizar modelos de reducción de dimensionalidad: si tenemos cientos (o miles, o incluso millones) de features, la reducción de dimensionalidad permite generar combinaciones de todas ellas, generando un dataset de variables virtuales mucho más compacto, más compatible con la mayoría de clasificadores y regresores. Los modelos de reducción de dimensionalidad son modelos de aprendizaje no supervisado (igual que los de clustering).
  3. Seleccionar variables, utilizando los propios resultados arrojados por los modelos. De esta forma podemos ir decidiendo si añadir o quitar más o menos variables, en función de si los resultados van mejorando o empeorando.

Realmente no sueles limitarte a tomar uno de estos tres caminos (¡sobre todo si es el primero!). Lo normal es usar una combinación de estas tres técnicas.

Nota: los modelos de reducción de dimensionalidad reducen drásticamente la explicabilidad del modelo. A veces esto no es relevante; pero otras veces sí lo es. Depende mucho del negocio y del valor que se quiere que aporte el machine learning. Si que un modelo sea opaco o difícil de entender no es problema, entonces bien. Pero en otros ámbitos (por ejemplo el bancario, donde si no ofreces un crédito a una persona debes decirle por qué) la explicabilidad y transparencia del modelo es esencial. En estos casos, la reducción de dimensionalidad no suele ser una opción viable, y el camino a elegir suele ser la selección de variables.

En este apartado vamos a explorar la reducción de dimensionalidad y la selección de variables.

8.A. Reducción de dimensionalidad

Si bien no hay una definición formal que englobe a los algoritmos de reducción de dimensionaliad (al menos que yo conozca), dichos algoritmos son fácilmente reconocibles: a partir de una serie de variables o features, genera otras completamente diferentes e independientes, intentanto extraer la relevancia y cualidades determinantes de las originales. El resultado es una menor cantidad de variables no originales (conocidas a veces como variables virtuales), pero muy significativas.

Estas variables virtuales generadas no tienen significado por sí mismas; son solo features generadas y que se asume que son significativas y útiles.

El modelo más famoso y utilizado en reducción de dimensionalidad es el PCA (Principal Component Analysis) y sus variantes, disponibles en sklearn.decomposition.

8.A.a. Principal Component Analysis

PCA o Principal Component Analysis transforma las variables o features que pasemos a dicho modelo en unas cuantas features virtuales (menos que las originales). Podemos pensar en PCA como en un pliegue de las variables en la dirección de la máxima varianza posible. De esta forma, las features virtuales generadas tendrán la mayor varianza posible; lo cual suele dar mayor poder discriminador al análisis que hagamos a posteriori. Supongamos que tenemos el siguiente dataset, con dos features $x_1$ y $x_2$:

pca_0

Y queremos obtener una única feature virtual a partir de esas dos. PCA va a intentar plegar los puntos en la dirección en la cual la varianza de la nueva feature resultante es la máxima posible.

Supongamos que realizamos el siguiente pliegue:

pca_1

La línea azul muestra el "nuevo eje de coordenadas" de la nueva variable generada. Como solo queríamos una, solo habrá un eje. No obstante, podemos ver que la distribución de esta nueva variable generada (los puntos azules) no tiene demasiada varianza; todos los puntos están muy cerca unos de otros. Esto es porque no hemos plegado en la dirección de la máxima varianza de las features originales.

Si de verdad utilizamos PCA, el resultado sería algo así:

pca_2

Ahora sí hemos plegado las dos features originales en la dirección de la máxima varianza posible. Los puntos azules tienen ahora mucha más varianza. Intuitivamente podemos ver que esta nueva feature virtual generada captura la mayoría de la varianza de las features originales.

Y esto es exactamente lo que hace PCA: intentar generar variables virtuales que aporten toda la información posible sobre las variables originales, pero en una representación más compacta (en este caso pasando de tener dos variables a solo una). Esto se llama utilizar la top principal component del dataset original.

En un caso donde tuviéramos originalmente 300 variables y quisiéramos quedarnos con 5, usaríamos PCA para quedarnos con las top 5 principal components del dataset.

Vamos a utilizar sklearn.decomposition.PCA sobre el dataset de clasificación que habíamos generado en el apartado de regresión. Era este (ahora mismo las labels nos dan igual):

In [119]:
dataset_clasificacion.plot(kind="scatter",
                           x="hormona_a",
                           y="hormona_b",
                           edgecolors="black",
                           linewidths=1, 
                           s=30,
                           figsize=(8,6))
pass

Apliquemos PCA para obtener la top (1) principal component del dataset:

In [120]:
from sklearn.decomposition import PCA

pca = PCA(n_components=1) # La top 1 principal component
         

# PCA tiene los métodos .fit() para entrenar,
# .transform() para transformar el dataset y obtener
# nuestras top components, y .fit_transform() para 
# hacer ambas en un solo paso:
top_principal_component = pca.fit_transform(dataset_clasificacion[["hormona_a", "hormona_b"]])

top_principal_component es un array con la nueva variable virtual generada (la top principal component). Vamos a meterla en un DataFrame, junto con las dos features originales:

In [121]:
dataset_clasificacion_con_pca = dataset_clasificacion.copy() # Así podemos copiar un DataFrame

dataset_clasificacion_con_pca["top 1 principal component"] = top_principal_component

# Reordenamos columnas para que quede más mono:
dataset_clasificacion_con_pca = dataset_clasificacion_con_pca[["hormona_a", "hormona_b", "top 1 principal component", "label"]]

dataset_clasificacion_con_pca[:10]
Out[121]:
hormona_a hormona_b top 1 principal component label
0 4.574208 5.609926 -1.680521 1
1 3.333068 3.301306 0.920645 0
2 4.452082 1.453083 1.781268 1
3 3.145367 3.635088 0.756899 0
4 3.722375 3.588740 0.461008 0
5 4.856292 4.993343 -1.340650 1
6 3.301592 4.371934 0.065430 1
7 5.804936 5.768588 -2.521728 1
8 5.694449 2.112604 0.524729 1
9 1.139096 5.987233 -0.001695 0

Vamos a ver la varianza:

In [122]:
dataset_clasificacion_con_pca.var() # Método .var() del DataFrame
Out[122]:
hormona_a                    1.489491
hormona_b                    1.607056
top 1 principal component    1.725830
label                        0.250417
dtype: float64

Podemos ver que la top principal component intenta capturar toda la varianza posible de hormona_a y hormona_b, dando como resultado una variable con más varianza que las dos originales. Ahora podríamos utilizar esta principal component como variable, en vez de las dos originales; de esta forma, tendríamos un problema con una dimensionalidad menor (1-D en vez de 2-D).

PCA también se utiliza mucho para visualización: imagínate que tenemos un dataset con 500 features, y queremos dibujar un scatterplot... Pues es imposible: solo podemos representar hasta 3 dimensiones de forma fácilmente entendible. Pues bien, podríamos utilizar PCA para obtener las top 3 principal components del dataset, y dibujar el scatterplot de los datos en función de estas 3 dimensiones creadas.

8.B. Selección de variables

Los modelos de reducción de dimensionalidad como PCA no deben ser utilizados a la ligera. Como hemos comentado antes, las variables generadas son virtuales, y carecen de significado como tal. En determinados ámbitos esto es aceptable (por ejemplo en neurociencia, donde las variables originales son zonas del cerebro; y hay miles de estas). Pero en otros ámbitos más estándar (marketing, banca y seguros, sector media...) normalmente queremos modelar con variables reales.

En estos casos, si queremos trabajar con menos variables, suele tocar seleccionar.

Para este ejemplo vamos a utilizar el dataset de los precios de las casas de Boston:

In [123]:
boston_training_set[:5]
Out[123]:
CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX PTRATIO B LSTAT Precio
384 20.08490 0.0 18.10 0.0 0.700 4.368 91.2 1.4395 24.0 666.0 20.2 285.83 30.63 8.8
205 0.13642 0.0 10.59 0.0 0.489 5.891 22.3 3.9454 4.0 277.0 18.6 396.90 10.87 22.6
313 0.26938 0.0 9.90 0.0 0.544 6.266 82.8 3.2628 4.0 304.0 18.4 393.39 7.90 21.6
50 0.08873 21.0 5.64 0.0 0.439 5.963 45.7 6.8147 4.0 243.0 16.8 395.56 13.45 19.7
90 0.04684 0.0 3.41 0.0 0.489 6.417 66.1 3.0923 2.0 270.0 17.8 392.18 8.81 22.6

Donde tenemos 13 features más nuestro target.

sklearn.feature_selection incluye varias técnicas para hacer selección de variables. Vamos a probar algunas.

Pero primero vamos a usar el poder de los DataFrames jerárquicos para hacer nuestras vida más sencilla. Vamos a generar una jerarquía de columnas que nos permita separar las features del target fácilmente:

In [124]:
columnas_jerarquicas = pd.MultiIndex.from_tuples([("features", "CRIM"),
                                                  ("features", "ZN"),
                                                  ("features","INDUS"),
                                                  ("features","CHAS"),
                                                  ("features","NOX"),
                                                  ("features","RM"),
                                                  ("features", "AGE"),
                                                  ("features", "DIS"),
                                                  ("features", "RAD"),
                                                  ("features", "TAX"),
                                                  ("features", "PTRATIO"),
                                                  ("features", "B"),
                                                  ("features", "LSTAT"),
                                                  ("target", "Precio")])
# Y las aplicamos:
boston_training_set.columns = columnas_jerarquicas

boston_training_set[:5]
Out[124]:
features target
CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX PTRATIO B LSTAT Precio
384 20.08490 0.0 18.10 0.0 0.700 4.368 91.2 1.4395 24.0 666.0 20.2 285.83 30.63 8.8
205 0.13642 0.0 10.59 0.0 0.489 5.891 22.3 3.9454 4.0 277.0 18.6 396.90 10.87 22.6
313 0.26938 0.0 9.90 0.0 0.544 6.266 82.8 3.2628 4.0 304.0 18.4 393.39 7.90 21.6
50 0.08873 21.0 5.64 0.0 0.439 5.963 45.7 6.8147 4.0 243.0 16.8 395.56 13.45 19.7
90 0.04684 0.0 3.41 0.0 0.489 6.417 66.1 3.0923 2.0 270.0 17.8 392.18 8.81 22.6

Mejor. Vamos a ello.

8.B.a. Filtro de variables por varianza

sklearn.feature_selection.VarianceThreshold nos elimina automáticamente las variables con baja varianza. La intuición es sencilla: variables que no varían mucho a lo largo del dataset, son poco informativas:

In [125]:
from sklearn.feature_selection import VarianceThreshold

filtro_var = VarianceThreshold(threshold=5.0) # Todas las features con menor varianza que
                                              # este número, serán filtradas.

# Lo "entrenamos" (no es un entrenamiento como tal,
# pero configura el selector de variables para nuestro
# dataset):
filtro_var.fit(boston_training_set["features"])

# Y transformamos nuestro dataset, efectivamente descartando
# variables:
boston_var_filtered = filtro_var.transform(boston_training_set["features"])

# La primera fila tras la transformación:
boston_var_filtered[0]
Out[125]:
array([  20.0849,    0.    ,   18.1   ,   91.2   ,   24.    ,  666.    ,
        285.83  ,   30.63  ])

Ha dejado nuestras 13 features en 8... Pero resulta complicado (o tedioso) saber cuáles ha filtrado, puesto que hemos perdido los nombres de las columnas.

Esta función permite recuperarlas:

In [126]:
def recuperar_columnas(dataframe_original, input_array):
    resultado = pd.DataFrame()
    input_dataframe = pd.DataFrame(input_array)
    input_dataframe.columns = [str(col) for col in input_dataframe.columns]
    for col_orig in dataframe_original.columns:
        for col_reduc in input_dataframe.columns:
            if np.all(dataframe_original[col_orig].values==input_dataframe[col_reduc].values):
                resultado[col_orig] = dataframe_original[col_orig]
            else:
                continue
    # Si es un DataFrame jerárquico:
    if type(resultado.columns[0])==tuple:
        multiIndex = pd.MultiIndex.from_tuples(resultado.columns)
        resultado.columns = multiIndex
    return resultado

La aplicamos:

In [127]:
boston_var_filtered = recuperar_columnas(boston_training_set, boston_var_filtered)

# Y volvemos a añadir el target:
boston_var_filtered["target","Precio"] = boston_training_set["target","Precio"]
boston_var_filtered[:5]
Out[127]:
features target
CRIM ZN INDUS AGE RAD TAX B LSTAT Precio
384 20.08490 0.0 18.10 91.2 24.0 666.0 285.83 30.63 8.8
205 0.13642 0.0 10.59 22.3 4.0 277.0 396.90 10.87 22.6
313 0.26938 0.0 9.90 82.8 4.0 304.0 393.39 7.90 21.6
50 0.08873 21.0 5.64 45.7 4.0 243.0 395.56 13.45 19.7
90 0.04684 0.0 3.41 66.1 2.0 270.0 392.18 8.81 22.6

El método de filtrar variables por varianza está bien, pero se basa bastante el el criterio personal de: ¿qué es una varianza alta?

8.B.b. Selección de variables univariante

sklearn presenta varias funciones para seleccionar variables en función de pruebas realizadas a las variables una por una (es decir, sin medir interacciones de ningún modo entre ellas).

Uno de ellos es sklearn.feature_selection.SelectKBest, que selecciona $K$ variables (nosotros elegimos $K$) a partir de una función de scoring que devuelva p-valores:

In [128]:
from sklearn.feature_selection import SelectKBest
from sklearn.feature_selection import f_regression

# Generamos la instancia de SelectKBest
kb = SelectKBest(score_func = f_regression, # Pasamos únicamente el método, sin argumentos
                 k=4 # Queremos quedarnos con las 4 variables mejores
                )

Y ahora utilizamos .fit_transform por ejemplo, para realizar la selección:

In [129]:
boston_KBest = kb.fit_transform(X=boston_training_set["features"],
                                y=boston_training_set["target","Precio"])

# Y lo pasamos a DataFrame:
boston_KBest = recuperar_columnas(boston_training_set, boston_KBest)

# Y re-añadimos la columna de target:
boston_KBest["target","Precio"] = boston_training_set["target","Precio"]

boston_KBest[:5]
Out[129]:
features target
INDUS RM PTRATIO LSTAT Precio
384 18.10 4.368 20.2 30.63 8.8
205 10.59 5.891 18.6 10.87 22.6
313 9.90 6.266 18.4 7.90 21.6
50 5.64 5.963 16.8 13.45 19.7
90 3.41 6.417 17.8 8.81 22.6

8.B.c. Uso de árboles u otros clasificadores y regresores para seleccionar variables

Como ya hemos visto, los árboles de decisión van realizando splits en función de las distintas variables para clasificar. Dichos splits comienzan por la variable más discriminante, y luego la siguiente, y luego la siguiente... Hasta que el árbol clasifica correctamente. Teniendo el orden de las variables por el que el árbol va haciendo los splits, tenemos el orden de importancia de las variables (desde el punto de vista de la opinión del árbol, que evidentemente no tiene por qué ser perfecto).

Los árboles no son los únicos modelos capaces de hacer esto. De hecho, cualquier modelo que tras ser entrenado con .fit() tenga los atributos .coef_ o .feature_importances_ puede ser utilizado para seleccionar variables.

La selección de variables acorde a un clasificador/regresor se puede hacer con sklearn.feature_selection.SelectFromModel:

In [130]:
from sklearn.feature_selection import SelectFromModel

# Creamos un arbol, que como el problema es de regresión...
arbol = DecisionTreeRegressor(criterion="mse", # Para que la métrica de criterio sea MSE
                                 max_depth=6,
                                 min_samples_split=5,
                                 min_samples_leaf=5,
                                 random_state=4
                                )

# Y se lo pasamos a una instancia de SelectFromModel:
sfm = SelectFromModel(estimator=arbol)

# Y usamos .fit_transform() para seleccionar:
boston_selectedFromModel = sfm.fit_transform(X=boston_training_set["features"],
                                             y=boston_training_set["target","Precio"])

# Lo pasamos a DataFrame y añadimos columna de target:
boston_selectedFromModel = recuperar_columnas(boston_training_set, boston_selectedFromModel)

# Y re-añadimos la columna de target:
boston_selectedFromModel["target","Precio"] = boston_training_set["target","Precio"]

boston_selectedFromModel[:5]
Out[130]:
features target
RM LSTAT Precio
384 4.368 30.63 8.8
205 5.891 10.87 22.6
313 6.266 7.90 21.6
50 5.963 13.45 19.7
90 6.417 8.81 22.6

Parece que el árbol ha decidido quedarse únicamente esas variables...

Ahora podríamos utilizar cualquier algoritmo (¡no necesariamente otro árbol!) para hacer nuestro modelo.

Recuerdo: el árbol construido y pasado a SelectFromModel solo se ha utilizado para seleccionar variables según su criterio.

8.B.d. Recursive Feature Elimination

La Recursive Feature Elimination es una de las técnicas de backward elimination para selección de variables. Comienzas con todas ellas, y pasas un clasificador (o regresor). Aquella variable que el modelo considera menos significativa es eliminada. Y así continuamente hasta que:

  1. Nos quedamos con el número de variables deseado.
  2. El modelo "decide" que no debe quitar ninguna más.

Para obtener el primer resultado (quedarnos con el número de features que queremos) utilizamos sklearn.feature_selection.RFE. Para obtener el segundo resultado se utiliza sklearn.feature_selection.RFECV: un Recursive Feature Elimination que además realiza validación cruzada para decidir cuántas variables debe dejar (más adelante explicaremos lo que es la validación cruzada).

Probemos los dos. Primero RFE:

In [131]:
from sklearn.feature_selection import RFE, RFECV

# Como modelo para realizar las RFEs vamos a usar un
# Support Vector Regressor:
mi_svr = SVR(kernel="linear")

# Y ahora instanciamos una RFE:
rfe = RFE(estimator = mi_svr,
          n_features_to_select=7, # Quiero que se quede con 7
          step=2 # Que vaya quitando variables de 2 en 2
         )

# Y utilizamos .fit_transform() para ejecutar el RFE y que
# seleccione:
boston_selected_rfe = rfe.fit_transform(X=boston_training_set["features"],
                                             y=boston_training_set["target","Precio"])

# Y a DataFrame(y metemos de nuevo el target):
boston_selected_rfe = recuperar_columnas(boston_training_set, boston_selected_rfe)
boston_selected_rfe["target","Precio"] = boston_training_set["target","Precio"]

boston_selected_rfe[:5]
Out[131]:
features target
INDUS CHAS NOX RM DIS PTRATIO LSTAT Precio
384 18.10 0.0 0.700 4.368 1.4395 20.2 30.63 8.8
205 10.59 0.0 0.489 5.891 3.9454 18.6 10.87 22.6
313 9.90 0.0 0.544 6.266 3.2628 18.4 7.90 21.6
50 5.64 0.0 0.439 5.963 6.8147 16.8 13.45 19.7
90 3.41 0.0 0.489 6.417 3.0923 17.8 8.81 22.6

Pidiéndole que nos dejara 7, el SVR ha decidido que deben ser esas 7. Pero... ¿y si no le decimos cuántas queremos, sino que le dejemos que él mismo decida? Para eso, utilizamos RFECV:

In [132]:
# Seguimos con el mismo SVR de antes, así que nada:
# creamos ya la instancia de RFECV:

rfecv = RFECV(estimator = mi_svr,
              step=1, # Que vaya quitando de una en una
              cv=10, # Número de "folds" de validación cruzada (lo explicamos luego)
              scoring="neg_mean_squared_error" # Métrica de rendimiento. Un string de las de esta lista:
                                               # http://goo.gl/bPGYAo
                                               # o bien una función de las disponibles en sklearn.metrics 
                                               # (de las métricas de regresión en este caso, claro).
              )

# Y utilizamos .fit_transform() para ejecutar el RFECV y que
# seleccione:
boston_selected_rfecv = rfecv.fit_transform(X=boston_training_set["features"],
                                             y=boston_training_set["target","Precio"])

# Y a DataFrame(y metemos de nuevo el target):
boston_selected_rfecv = recuperar_columnas(boston_training_set, boston_selected_rfecv)
boston_selected_rfecv["target","Precio"] = boston_training_set["target","Precio"]

boston_selected_rfecv[:5]
Out[132]:
features target
INDUS CHAS NOX RM DIS PTRATIO LSTAT Precio
384 18.10 0.0 0.700 4.368 1.4395 20.2 30.63 8.8
205 10.59 0.0 0.489 5.891 3.9454 18.6 10.87 22.6
313 9.90 0.0 0.544 6.266 3.2628 18.4 7.90 21.6
50 5.64 0.0 0.439 5.963 6.8147 16.8 13.45 19.7
90 3.41 0.0 0.489 6.417 3.0923 17.8 8.81 22.6

Y decidió quedarse con esas.

Recomiendo prestar siempre muhcísima atención a la selección de variables. RFE y RFECV suelen funcionar muy bien. También es costumbre probarlos con diferentes clasificadores/regresores. Si todos ellos coinciden en determinada selección, es probablemente porque es la más adecuada.

Existen muchas técnicas más para la selección de variables. Una de ellas es la eliminación bidireccional o stepwise; que es similar a RFE, con la diferencia de que va eliminando pero también añadiendo variables, hasta probar todas las combinaciones posibles. Es computacionalmente muy cara (entrena el modelo para todas las posibles combinaciones de variables), pero resulta sencilla de implementar con un par de bucles for y los clasificadores y métricas de rendimiento que nos ofrece sklearn.

Dicho esto, pasemos a la validación y selección de modelos.

9. Validación y selección de modelos

Supongamos que ya hemos seleccionado variables para lanzar modelos. Conocemos cuáles hay, y sabemos más o menos cómo funcionan. Pero... ¿y todos los argumentos e hiperparámetros que toma cada modelo? ¿Cómo elijo eso?

La respuesta, tanto para seleccionar un modelo como sus hiperparámetros y configuraciones, es la validación.

La validación consiste en la aplicación de métricas para saber qué modelo (o variante de modelo) funciona mejor sobre datos que no ha visto antes. Parece que hace lo mismo que la separación en training y testing, ¿no?

Sí y no. Existen varios enfoques para la validación y testing de modelos. El más utilizado es el siguiente. Se parece a una especie de juego:

  1. Dividimos el dataset en dos conjuntos: el de training y testing. El conjunto de testing lo guardamos en una caja fuerte, y no vamos a volver a verlo ni utilizarlo hasta el final del proceso.
  2. Quedándonos con el conjunto de training, tenemos que hacer lo que podamos para elegir un modelo. Solo podemos elegir uno: tanto el modelo en sí como sus hiperparámetros y configuraciones.
  3. Cuando ya tengamos el modelo elegido (por ejemplo: regresión logística con regularización L2 y alpha de 0.1, con fit_intercept=True), lo entrenamos con TODO nuestro training set.
  4. Ahora ya, podemos abrir nuestra caja fuerte con el conjunto de testing. Los resultados de precisión obtenidos en este conjunto es lo que podemos esperar de nuestro modelo para el futuro.

Intuirás que lo complicado es el paso 2. Ahí es donde entra en juego la validación.

En validación, cogemos parte de nuestro training set y lo utilizamos para medir el rendimiento de nuestros modelos. Nos quedaremos con el que mejor rendimiento dé. De forma que lo que hemos ido haciendo en este Notebook ha sido una forma muy informal de validación; a pesar de haberla hecho con un conjunto de datos que hemos llamado testing (lo llamé testing para no meter al principio del todo el tema de la validación)

Repito: para ser léxicamente correctos: el conjunto de testing no debería verse ni utilizarse hasta que hayamos seleccionado nuestro modelo ganador mediante validación; y este conjunto de validación lo cogemos de una parte del conjunto de training.

9.A. Training, validación, testing

Hay distintas técnicas de validación. La más básica consiste en separar una parte del conjunto de training y llamarla conjunto de validación, de forma que:

  1. Entrenamos modelos con el conjunto de training.
  2. Comparamos, y seleccionamos el modelo que mejor rendimiento nos da en el conjunto de validación.
  3. Utilizamos training + validación para entrenar de nuevo el modelo ganador, sobre los dos conjuntos de datos a la vez.
  4. Medimos el rendimiento en el conjunto de testing; lo ponemos en una PPT y se lo entregamos a los inversores.

Para hacerlo, podemos utilizar sklearn.cross_validation.train_test_split. Lo hacemos una vez para separar en training y testing, y luego otra vez con el conjunto de training para separar en training de verdad y validación.

Lioso, lo sé. Este dibujito lo muestra de forma más sencilla: training_val_test

Siempre surge la pregunta: ¿qué cantidad de datos dejar para cada conjunto? Existen varios compromisos:

  1. Si dejamos pocos datos para el test set y validation set, igual no son representativos del todo...
  2. Pero si cogemos muchos, los quitamos de training y nuestro modelo entrenará con menos datos...

Un criterio bastante utilizado es 60-20-20 (60% entrenamiento, 20% validación y 20% test). No funciona mal... Depende sobre todo de la cantidad de datos que tengamos. En general: cuantos menos datos, peor funciona este método (y todos; pero éste en especial).

Por suerte existe una forma alternativa, que permite que el conjunto de validación sea pequeño y grande a la vez... Y ahí entra en juego la validación cruzada.

9.B. Cross Validation

Cross Validation o validación cruzada son un conjunto de técnicas que hacen separaciones dinámicas y automáticas de conjuntos de training y validación.

Hay varias variantes de validación cruzada. Vamos a centrarnos en la más utilizada: K-Fold Cross Validation.

K-Fold Cross Validation consiste en que por cada modelo, no entrenamos una vez, sino un número $K$ de veces. Pongamos un valor de $K$ de 10. Pues bien, entrenaremos nuestro modelo 10 veces, cada vez dejando un 10% de los datos para validación; y cada vez siendo este 10% distinto. Una vez entrenado el modelo las 10 veces, sacamos el rendimiento del modelo haciendo la media de los resultados obtenidos ($MSE$, porcentaje bien clasificado, etc) en cada una de estas 10 veces. De esta forma, a nivel efectivo hemos entrenado con el 100% de los datos de training; y también hemos validado con el 100% de los datos de training (de nuevo: no tocamos el conjunto de testing).

Esta imagen (por cortesía de Neerav Basant) lo ilustra muy bien:

10-fold-cv

Podemos elegir el valor de $K$ que queramos. Lo normal es un valor entre 3 y 10. Cuanto mayor sea la $K$, más modelos serán entrenados, y más lento será el proceso (pero más robusto en teoría).

Exploremos ahora las distintas formas de hacer validación cruzada con sklearn:

9.B.a. Validación cruzada empotrada en la especificación de los modelos

Ciertos modelos de sklearn (principalmente los de sklearn.linear_model) tienen una "variante" CV, que nos permite elegir los mejores hiperparámetros y configuraciones de una forma sencilla.

Por ejemplo: la regresión Ridge tiene su variante RidgeCV que nos permite pasarle a la especificación del modelo una serie de hiperparámetros $\alpha$ a explorar (el hiperparámetro de shrinkage de la regularización $L2$). Sklearn entrenará cada modelo $K$ veces haciendo validación cruzada y nos devolverá el mejor modelo entrenado. Utilicémoslo con nuestro dataset de training de regresión (de años trabajados ~ salario anual):

In [133]:
from sklearn.linear_model import RidgeCV

ridge_cv = RidgeCV(alphas=[0.05, 0.1, 0.5, 1.0, 5.0, 10.0, 20.0], # Queremos probar con estos valores
                   fit_intercept=True,
                   cv=10 # 10-Fold Cross Validation
                   )

# Lo entrenamos. En el fondo, estamos entrenando
# 10 folds * 8 valores de alpha distintos = 80 modelos;
# así que tarda un poquito:
ridge_cv.fit(X=dataset_regresion_training[["años trabajados"]],
             y=dataset_regresion_training["salario bruto"])
Out[133]:
RidgeCV(alphas=[0.05, 0.1, 0.5, 1.0, 5.0, 10.0, 20.0], cv=10,
    fit_intercept=True, gcv_mode=None, normalize=False, scoring=None,
    store_cv_values=False)

Modelos entrenados, y sklearn se ha "quedado" automáticamente con el mejor. El mejor de ellos ha sido el que tiene el hiperparámetro $\alpha$:

In [134]:
ridge_cv.alpha_
Out[134]:
20.0

Ese es el hiperparámetro de los elegidos que minimiza el $MSE$ en validación cruzada.

Los modelos que tienen una variante CV consiguen hacer la validación cruzada de una forma muy eficiente, pero son poco flexibles (normalmente te dejan elegir los valores a explorar de uno o dos hiperparámetros, pero no de todo).

Si queremos trabajar más la validación cruzada, podemos utilizar directamente sklearn.cross_validation.KFold, que genera un iterador sobre los datos (para que usemos un bucle for, vayamos entrenando con distintos valores de los hiperparámetros y argumentos que queramos, y vayamos acumulando los resultados en una lista/diccionario por ejemplo).

Pero aún mejor: GridSearchCV nos permite hacer lo mismo de una forma más elegante y menos tediosa.

9.B.b. GridSearchCV

Grid Search es una técnica que nos permite construir una gran cantidad de modelos, entrenarlos y validarlos (con validación cruzada) de forma sencilla y expresiva.

El Grid Search está compuesto por:

  • Un modelo (supervisado: clasificador o regresor)
  • Un grid de hiperparámetros, donde definimos para cada hiperparámetro los valores que queremos probar
  • Una técnica de validación cruzada
  • Una métrica de rendimiento

Al igual que en los modelos con validación cruzada empotrada, Grid Search entrena los modelos con validación cruzada y probando, en este caso, las distintas combinaciones de hiperparámetros; y midiendo el rendimiento de cada modelo. El resultado final es el modelo que mejor rendimiento ha ofrecido.

GridSearchCV es la clase de sklearn que nos permite hacer Grid Search.

Vamos a probar en nuestro dataset de clasificación (calvos/no-calvos). Supongamos que queremos probar con K-Nearest Neighbors: a priori no sabemos qué valor de $K$ será mejor, así que podemos construir una Grid Search. El grid de hiperparámetros lo construimos con un diccionario, donde cada clave es un string con el nombre exacto de cada hiperparámetro (en el caso de sklearn.neighbors.KNeighborsClassifier es "n_neighbors"). Si el modelo toma más hiperparámetros/argumentos de los cuales no especificamos valores a explorar, los dejará constantes (y serán los que vienen por defecto en el modelo):

In [135]:
from sklearn.model_selection import GridSearchCV

knn_gs = GridSearchCV(estimator = KNeighborsClassifier(), # Instanciamos un KNN
                      param_grid = {"n_neighbors":[1,3,5,7,9,11]}, # Los valores de n_neighbors a probar, en lista
                      scoring="accuracy", # La métrica de rendimiento a utilizar
                      cv=5, # 5-Fold Cross Validation
                      verbose=3 # Esto hace que imprima "cosas informativas" durante la ejecución
                     )

knn_gs es ahora una instancia de GridSearchCV, que podemos tratar como un modelo normal y corriente (con sus métodos .fit() y .predict(), por ejemplo). Si aplicamos el método .fit():

In [136]:
knn_gs.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
           y=dataset_clasificacion_training["label"])
pass
Fitting 5 folds for each of 6 candidates, totalling 30 fits
[CV] n_neighbors=1 ...................................................
[CV] .......................... n_neighbors=1, score=0.865979 -   0.0s
[CV] n_neighbors=1 ...................................................
[CV] .......................... n_neighbors=1, score=0.885417 -   0.0s
[CV] n_neighbors=1 ...................................................
[CV] .......................... n_neighbors=1, score=0.812500 -   0.0s
[CV] n_neighbors=1 ...................................................
[CV] .......................... n_neighbors=1, score=0.833333 -   0.0s
[CV] n_neighbors=1 ...................................................
[CV] .......................... n_neighbors=1, score=0.831579 -   0.0s
[CV] n_neighbors=3 ...................................................
[CV] .......................... n_neighbors=3, score=0.917526 -   0.0s
[CV] n_neighbors=3 ...................................................
[CV] .......................... n_neighbors=3, score=0.927083 -   0.0s
[CV] n_neighbors=3 ...................................................
[CV] .......................... n_neighbors=3, score=0.916667 -   0.0s
[CV] n_neighbors=3 ...................................................
[CV] .......................... n_neighbors=3, score=0.916667 -   0.0s
[CV] n_neighbors=3 ...................................................
[CV] .......................... n_neighbors=3, score=0.884211 -   0.0s
[CV] n_neighbors=5 ...................................................
[CV] .......................... n_neighbors=5, score=0.917526 -   0.0s
[CV] n_neighbors=5 ...................................................
[CV] .......................... n_neighbors=5, score=0.916667 -   0.0s
[CV] n_neighbors=5 ...................................................
[CV] .......................... n_neighbors=5, score=0.927083 -   0.0s
[CV] n_neighbors=5 ...................................................
[CV] .......................... n_neighbors=5, score=0.927083 -   0.0s
[CV] n_neighbors=5 ...................................................
[CV] .......................... n_neighbors=5, score=0.894737 -   0.0s
[CV] n_neighbors=7 ...................................................
[CV] .......................... n_neighbors=7, score=0.917526 -   0.0s
[CV] n_neighbors=7 ...................................................
[CV] .......................... n_neighbors=7, score=0.916667 -   0.0s
[CV] n_neighbors=7 ...................................................
[CV] .......................... n_neighbors=7, score=0.937500 -   0.0s
[CV] n_neighbors=7 ...................................................
[CV] .......................... n_neighbors=7, score=0.927083 -   0.0s
[CV] n_neighbors=7 ...................................................
[CV] .......................... n_neighbors=7, score=0.894737 -   0.0s
[CV] n_neighbors=9 ...................................................
[CV] .......................... n_neighbors=9, score=0.917526 -   0.0s
[CV] n_neighbors=9 ...................................................
[CV] .......................... n_neighbors=9, score=0.927083 -   0.0s
[CV] n_neighbors=9 ...................................................
[CV] .......................... n_neighbors=9, score=0.937500 -   0.0s
[CV] n_neighbors=9 ...................................................
[CV] .......................... n_neighbors=9, score=0.937500 -   0.0s
[CV] n_neighbors=9 ...................................................
[CV] .......................... n_neighbors=9, score=0.905263 -   0.0s
[CV] n_neighbors=11 ..................................................
[CV] ......................... n_neighbors=11, score=0.917526 -   0.0s
[CV] n_neighbors=11 ..................................................
[Parallel(n_jobs=1)]: Done   1 out of   1 | elapsed:    0.0s remaining:    0.0s
[Parallel(n_jobs=1)]: Done   2 out of   2 | elapsed:    0.0s remaining:    0.0s
[CV] ......................... n_neighbors=11, score=0.916667 -   0.0s
[CV] n_neighbors=11 ..................................................
[CV] ......................... n_neighbors=11, score=0.937500 -   0.0s
[CV] n_neighbors=11 ..................................................
[CV] ......................... n_neighbors=11, score=0.916667 -   0.0s
[CV] n_neighbors=11 ..................................................
[CV] ......................... n_neighbors=11, score=0.905263 -   0.0s
[Parallel(n_jobs=1)]: Done  30 out of  30 | elapsed:    0.2s finished

Nos ha impreso cómo ha ido entrenando y validando cada modelo con validación cruzada, el valor de n_neighbors que ha probado cada vez, el rendimiento obtenido y el tiempo que ha tardado en entrenar.

Podemos preguntarle a knn_gs cuál ha sido el mejor modelo en validación con el atributo .best_estimator_:

In [137]:
knn_gs.best_estimator_
Out[137]:
KNeighborsClassifier(algorithm='auto', leaf_size=30, metric='minkowski',
           metric_params=None, n_jobs=1, n_neighbors=9, p=2,
           weights='uniform')

Una vez entrenado, validado y elegido con el método .fit(), knn_gs utilizará el modelo ganador si aplicamos knn_gs.predict(), de forma que no es necesario que creemos otro KNeighborsClassifier con el número n_negihbors ganador para empezar a predecir con el mejor modelo. No obstante, lo recomiendo: utiliza GridSearchCV para conseguir la información de los hiperparámetros que mejor funcionan. Una vez que los conozcas, entrena el modelo de la forma convencional sin GridSearchCV; con los mejores hiperparámetros y argumentos.

Hagamos otro ejemplo de GridSearchCV, esta vez con árboles regresores sobre nuestro dataset de regresión (años trabajados ~ salario anual). Lo bueno es que los árboles toman muchos argumentos e hiperparámetros distintos; y GridSearchCV entrena todas las combinaciones posibles de todos ellos (que puede ser lento, pero es eficaz):

In [138]:
# Creamos el Grid Search de árboles:
arbol_gs = GridSearchCV(estimator=DecisionTreeRegressor(),
                        param_grid={"max_depth": [2,3,4,5,6],
                                    "min_samples_split": [2,5,10,15],
                                    "min_samples_leaf": [2,5,10,15]},
                        scoring="neg_mean_squared_error",
                        cv=5,
                        verbose=1 # Para que imprima menos información que antes
                       )

# Entrenamos, validamos y elegimos el mejor árbol:
arbol_gs.fit(X=dataset_regresion_training[["años trabajados"]],
             y=dataset_regresion_training["salario bruto"]) 
pass
Fitting 5 folds for each of 80 candidates, totalling 400 fits
[Parallel(n_jobs=1)]: Done 400 out of 400 | elapsed:    2.0s finished

Vemos que a medida que metemos más hiperparámetros y más valores a explorar para cada uno de ellos, el número de modelos entrenados (y validados) aumenta exponencialmente. Veamos cuál es nuestro árbol ganador:

In [139]:
arbol_gs.best_estimator_
Out[139]:
DecisionTreeRegressor(criterion='mse', max_depth=5, max_features=None,
           max_leaf_nodes=None, min_impurity_split=1e-07,
           min_samples_leaf=2, min_samples_split=2,
           min_weight_fraction_leaf=0.0, presort=False, random_state=None,
           splitter='best')

GridSearchCV nos ayda muchísimo a elegir los mejores hiperparámetros y argumentos para un modelo. Pero... ¿Y cuando queremos elegir entre varios modelos distintos (por ejemplo, entre el mejor árbol y el mejor SVM)?

Sencillo: resulta que cualquier GridSearchCV generado nos devuelve también el rendimiento del mejor modelo a través del atributo .best_score_:

In [140]:
arbol_gs.best_score_
Out[140]:
-29566703.775282249

Que en este caso es el $MSE$ (puesto que así lo definimos en arbol_gs al generarlo). Que no te despiste el signo menos que tiene delante el $MSE$: es por razones internas de sklearn. Básicamente, en la versión 0.18 (la que estamos utilizando aquí) se calcula el $-MSE$ (en negativo). Pero es exactamente lo mismo; solo que en vez de minimizar el $MSE$, se maximiza el $-MSE$, que es completamente equivalente.

El caso es que podemos:

  • Construir varios GridSearchCV: uno para cada modelo. Esto nos permitirá elegir los mejores hiperparámetros para cada modelo.
  • Meter todos los GridSearchCV en una lista de Python normal y corriente.
  • Entrenar todos ellos con .fit(), y obtener los .best_score_ de cada modelo ganador.
  • Elegir el modelo con mejor .best_score_

Vamos a ponerlo en práctica con nuestro problema de clasificación de la alopecia. Vamos a probar los siguientes modelos:

  • Regresión logística
  • K-NN
  • SVM
  • Random Forest

Buscando los mejores hiperparámetros para cada uno, y luego quedándonos con el mejor modelo. Manos a la obra:

In [141]:
# Definimos todos los modelos con su GridSearchCV:

reg_log_gs = GridSearchCV(estimator=LogisticRegression(),
                          param_grid={"penalty": ["l1","l2"], # Probamos con regularizaciones L1 y L2
                                      "C": [0.5, 1.0, 1.5]}, # Y varios valores de la regularización C
                          scoring="roc_auc", # En vez de accuracy, vamos a usar el área bajo la curva ROC
                          cv=5)

knn_gs = GridSearchCV(estimator=KNeighborsClassifier(),
                      param_grid={"n_neighbors": [1,3,5,7,9,11]},
                      scoring="roc_auc", # Lo suyo es usar la misma métrica en todos los modelos que probemos
                      cv=5)

svm_gs = GridSearchCV(estimator=SVC(),
                      param_grid={"C": [0.5, 1.0, 1.5],
                                  "kernel": ["linear","poly","rbf"],
                                  "degree": [2,3,4]},
                      scoring="roc_auc",
                      cv=5)

random_f_gs = GridSearchCV(estimator=RandomForestClassifier(),
                           param_grid={"n_estimators": [5,10,20,50],
                                       "max_depth": [2,3,4,5],
                                       "min_samples_split": [3,4,5,10],
                                       "min_samples_leaf": [3,4,5,10]},
                           scoring="roc_auc",
                           cv=5)

# Los metemos todos en una lista:
bag_de_modelos = [reg_log_gs, knn_gs, svm_gs, random_f_gs]

# Los entrenamos todos, validando y eligiendo el mejor de cada categoría
# (la mejor regresión logística, el mejor K-NN, el mejor SVM y el mejor Random Forest):
for modelo in bag_de_modelos:
    modelo.fit(X=dataset_clasificacion_training[["hormona_a", "hormona_b"]],
               y=dataset_clasificacion_training["label"])
    
# Creamos otra lista para meter los .best_score_:
lista_de_scores = []

# Y los metemos:
for modelo in bag_de_modelos:
    lista_de_scores.append(modelo.best_score_)

Y los resultados son (recordamos que las listas de Python mantienen el orden de los elementos; así que el primer score es de la regresión logística, el segundo de K-NN, etc):

In [142]:
lista_de_scores
Out[142]:
[0.96781277892362616,
 0.96604229236020656,
 0.97501366348844487,
 0.96733889145197083]

Ha ganado el tercer modelo, es decir, el SVM (SVC): es el que mayor área bajo la curva ROC ha dado. Vamos a ver con qué argumentos/hiperparámetros ha ganado:

In [143]:
bag_de_modelos[2].best_estimator_
Out[143]:
SVC(C=1.0, cache_size=200, class_weight=None, coef0=0.0,
  decision_function_shape=None, degree=4, gamma='auto', kernel='poly',
  max_iter=-1, probability=False, random_state=None, shrinking=True,
  tol=0.001, verbose=False)

Y de los generados, este es el modelo que querríamos utilizar para hacer nuestras predicciones a futuro.

La validación, así como la comparación y selección de modelos es la esencia del aprendizaje automático. Siempre buscamos un modelo capaz de conseguir buenos resultados con datos con los que no ha entrenado; y la validación cruzada junto con la búsqueda de hiperparámetros nos permite encontrar el modelo con el menor $E_{out}$ posible.

Por último, vamos a ver un elemento que cada vez está cobrando más importancia en la ciencia de datos: las pipelines.

10. Pipelines

Ya hemos visto las principales ramas de sklearn. Solo queda por comentar el tema de las pipelines.

Una pipeline no es más que una especie de diagrama que encadena distintos procesos. Normalmente, tenemos que seguir un determinado orden a la hora de utilizar aprendizaje automático: preprocesar y transformar los datos, seleccionar variables, entrenar modelos, validarlos y seleccionarlos. Las pipelines intentan simplificarnos la tarea: encadenan varias de estas fases, y generan un modelo portable, fácil de entrenar y entender.

Las pipelines son utilizadas en muchas bibliotecas de machine learning. Si alguna vez has usado herramientas visuales de aprendizaje automático y minería de datos (como SAS Enterprise Miner, SPSS Modeller, RapidMiner o KNIME), seguro habrás "arrastrado" distintos transformadores y modelos, los habrás configurado, y luego los habrás conectado entre sí. Pues una pipeline en sklearn es exactamente lo mismo.

Sklearn implementa sklearn.pipeline.Pipeline para construir pipelines.

Vamos a utilizar el dataset del precio de las casas de Boston. Querremos hacer lo siguiente, y por el siguiente orden:

  1. Seleccionar variables
  2. Puesto que las variables son de variopintos ámbitos y magnitudes, probablemente sea buena idea estandarizarlas (o igual no...)
  3. Entrenar, validar y seleccionar modelos

Vamos a hacerlo todo con una Pipeline:

In [144]:
from sklearn.pipeline import Pipeline

# Lo primero será generar un RFECV para seleccionar variables,
# que utilizará un árbol regresor sencillo para ver qué variables
# son las que debemos seleccionar:
arbol_seleccionar = DecisionTreeRegressor()
rfecv = RFECV(estimator=arbol_seleccionar)

# Luego utilizaremos un StandardScaler para estandarizar
# las variables:
scaler = StandardScaler()

# Luego probaremos diferentes random forest regresores:
random_forest_r = RandomForestRegressor()

# Ahora generamos la Pipeline, donde pasaremos una lista
# de las "steps" o pasos a seguir, a modo de tuplas:
# donde el primer elemento es un string con el "nombre"
# de cada cosa, y el segundo elemento es cada cosa en sí:
mi_pipeline = Pipeline(steps=[("seleccionador",rfecv), ("escalador",scaler), ("forest",random_forest_r)])

mi_pipeline
Out[144]:
Pipeline(steps=[('seleccionador', RFECV(cv=None,
   estimator=DecisionTreeRegressor(criterion='mse', max_depth=None, max_features=None,
           max_leaf_nodes=None, min_impurity_split=1e-07,
           min_samples_leaf=1, min_samples_split=2,
           min_weight_fraction_leaf=0.0, presort=False, random_...timators=10, n_jobs=1, oob_score=False, random_state=None,
           verbose=0, warm_start=False))])

Hemos generado la pipeline. No obstante, vamos a querer hacer validación en todas las steps: en el modelo seleccionador de variables, en el escalador, y en el random forest. Para eso, podemos generar un diccionario que luego pasaremos a GridSearchCV, junto con la propia pipeline: utilizando los nombres personalizados que le hemos puesto a cada step, junto con la notación de poner doble underscore (__) y cada argumento/hiperparámetro, podemos definir los hiperparámetros que queremos probar de cada parte de la pipeline:

In [145]:
hiperparam_pipeline = {"seleccionador__cv":[5,10], # Si queremos explorar varios, los pasamos en lista
                       "seleccionador__scoring":["neg_mean_squared_error"], # Si no es así, lista con 1 elemento
                       "escalador__with_mean":[True, False], # Esta y la de abajo nos permiten explorar
                       "escalador__with_std":[True, False],  # el usar el escalador, o que "no haga nada";
                                                             # que es lo mismo que no ponerlo; de esta forma, 
                                                             # explorarmos si es mejor estandarizar o no
                                                             # en este caso. Si pones tanto with_mean como
                                                             # with_std en False en StandardScaler, es como
                                                             # si no estuviera, porque ni hace que la media
                                                             # de cada variable sea 0; ni que la desviación
                                                             # típica sea 1. También permite tocar la media 
                                                             # y no la varianza, o al revés. Es decir: todas
                                                             # las combinaciones.
                        "forest__n_estimators":[10,20,40], 
                        "forest__criterion":["mse"],
                        "forest__max_depth":[5,10,20],
                        "forest__min_samples_split":[2,3,5,10],
                        "forest__min_samples_leaf":[2,3,5]
                      }

Ahora podemos concebir mi_pipeline como "un gran modelo"; y por otro lado tenemos un diccionario con los hiperparámetros que queremos probar y validar en cada paso. Tal y como si fuera un modelo normal y corriente de sklearn, mi_pipeline tiene el método .fit(); que aplica .fit() en todos los elementos que tiene dentro.

Gracias a esto, podemos utilizar GridSearchCV para probar qué combinaciones son las mejores. Solo tenemos que pasar como modelo mi_pipeline, y como grid de hiperparámetros a explorar el diccionario hiperparametros_pipeline que hemos construido:

In [146]:
mi_pipeline_grid_search = GridSearchCV(estimator=mi_pipeline,
                                       param_grid=hiperparam_pipeline,
                                       scoring="neg_mean_squared_error",
                                       cv=10,
                                       n_jobs=4 # Para que use 4 CPUs a la vez; y entrene 
                                                # un modelo en cada una a la vez.
                                      )

# Y le damos caña. Esto va a ser largo:
mi_pipeline_grid_search.fit(X=boston_training_set["features"],
                            y=boston_training_set["target","Precio"])
Out[146]:
GridSearchCV(cv=10, error_score='raise',
       estimator=Pipeline(steps=[('seleccionador', RFECV(cv=None,
   estimator=DecisionTreeRegressor(criterion='mse', max_depth=None, max_features=None,
           max_leaf_nodes=None, min_impurity_split=1e-07,
           min_samples_leaf=1, min_samples_split=2,
           min_weight_fraction_leaf=0.0, presort=False, random_...timators=10, n_jobs=1, oob_score=False, random_state=None,
           verbose=0, warm_start=False))]),
       fit_params={}, iid=True, n_jobs=4,
       param_grid={'forest__criterion': ['mse'], 'seleccionador__cv': [5, 10], 'forest__min_samples_split': [2, 3, 5, 10], 'seleccionador__scoring': ['neg_mean_squared_error'], 'forest__min_samples_leaf': [2, 3, 5], 'escalador__with_std': [True, False], 'forest__n_estimators': [10, 20, 40], 'forest__max_depth': [5, 10, 20], 'escalador__with_mean': [True, False]},
       pre_dispatch='2*n_jobs', refit=True, return_train_score=True,
       scoring='neg_mean_squared_error', verbose=0)

Vamos a ver cuál es el modelo (o mejor dicho: pipeline) ganador con .best_estimator_:

In [147]:
pipeline_ganadora = mi_pipeline_grid_search.best_estimator_

pipeline_ganadora
Out[147]:
Pipeline(steps=[('seleccionador', RFECV(cv=10,
   estimator=DecisionTreeRegressor(criterion='mse', max_depth=None, max_features=None,
           max_leaf_nodes=None, min_impurity_split=1e-07,
           min_samples_leaf=1, min_samples_split=2,
           min_weight_fraction_leaf=0.0, presort=False, random_st...timators=20, n_jobs=1, oob_score=False, random_state=None,
           verbose=0, warm_start=False))])

Vamos a ver los steps de la pipeline ganadora en más detalle con el atributo .steps, que devuelve una lista con los 3 procesos (en este caso):

In [148]:
pipeline_ganadora.steps
Out[148]:
[('seleccionador', RFECV(cv=10,
     estimator=DecisionTreeRegressor(criterion='mse', max_depth=None, max_features=None,
             max_leaf_nodes=None, min_impurity_split=1e-07,
             min_samples_leaf=1, min_samples_split=2,
             min_weight_fraction_leaf=0.0, presort=False, random_state=None,
             splitter='best'),
     n_jobs=1, scoring='neg_mean_squared_error', step=1, verbose=0)),
 ('escalador', StandardScaler(copy=True, with_mean=True, with_std=False)),
 ('forest',
  RandomForestRegressor(bootstrap=True, criterion='mse', max_depth=20,
             max_features='auto', max_leaf_nodes=None,
             min_impurity_split=1e-07, min_samples_leaf=2,
             min_samples_split=10, min_weight_fraction_leaf=0.0,
             n_estimators=20, n_jobs=1, oob_score=False, random_state=None,
             verbose=0, warm_start=False))]

Y esa es la combinación ganadora. Ahora, lo ideal sería re-entrenar con los hiperparámetros que han resultado mejores, sobre todo el dataset de training (sin validación, porque ya la hemos hecho). Habiendo decidido que este es nuestro modelo (o pipeline) ganador, podríamos ver el rendimiento que da en el conjunto de testing; así tendríamos una idea de lo bien que debería predecir el modelo en el futuro.

Los modelos (así como los pre-procesamientos previos, selectores de variables...) deben ser re-entrenados con el tiempo. En el caso de que llegara una cantidad considerable de datos nuevos, lo ideal sería repetir el proceso desde 0. De esta forma "garantizamos" que nuestro trabajo sigue siendo útil a lo largo del tiempo.

11. Bibliografía y recursos adicionales

Dividiré esta sección en dos partes.

Por un lado, los libros que recomiendo relacionados con el aprendizaje automático y la minería de datos en general (no atados a Python, sino más a la matemática en la que se fundamenta cada modelo):

  • The Elements of Statistical Learning: EL libro de aprendizaje automático. Trata una gran variedad de modelos y técnicas, explicadas en gran detalle. Obligatorio. Punto.
  • Pattern Recognition and Machine Learning: también muy buen libro sobre aprendizaje automático. Eso sí: si eres muy frecuentista (y por tanto poco bayesiano) puedes acabar odiándolo. Recomendado a los amantes de la probabilidad y de Bayes, ya que intenta explicar todos los modelos a partir de ahí.
  • Data Mining: un libro bastante reciente (2015). No solo entra en los modelos, sino que presta bastante atención a las cosas que los rodean (como el feature engineering y las métricas de rendimiento de modelos y de distancia).
  • Bayesian Nonparametric Data Analysis: si quieres aprender más sobre cómo utilizar modelos bayesianos para el aprendizaje automático. Corto, pero intenso.

Y hay muchísimos más.

Luego ya dentro de Python hay buenos libros (como Building Machine Learning Systems with Python y Python Machine Learning) sobre machine learning con sklearn, además de otros más avanzados como Advanced Machine Learning with Python (sobre modelos menos convencionales y más avanzados) y Large Scale Machine Learning with Python (sobre machine learning a gran escala con Python) que son muy nuevos. Además, existen muchas bibliotecas más allá de sklearn útiles para el mundo del data science y machine learning:

  • Bibliotecas de redes neuronales y Deep Learning: las comentamos anteriormente. Destaco Theano, TensorFlow y Keras.
  • Bibliotecas para procesamiento de lenguaje natural: NLTK es una de las bibliotecas más famosas (si no la más) para hacer procesamiento de lenguaje natural. Ofrece muchas funciones y utilidades para el prepocesado de texto no estructurado (tokenizers, stemmers, computadores de n-grams, lemmatizers, part-of-speech taggers...) así como modelos de machine learning para utilizar una vez que el texto ha sido "traducido" a datos que los modelos pueden utilizar. No obstante, es frecuente utilizar NLTK para preprocesar el texto y utilizar luego sklearn para los modelos en sí (porque son más potentes). Genism es otra biblioteca muy buena; esencial para topic modelling y word embeddings (que incluye modelos como Word2vec). Los mejores libros para hacer procesamiento de lenguaje natural con Python son Natural Language Processing with Python (segunda edición en camino) y Python 3 Text Processing with NLTK 3 Cookbook.
  • Bibliotecas para visión computacional y procesamiento de imágenes: OpenCV es la biblioteca más famosa para visión computacional: escrita en C++, incluye API para Python y hay buenos libros sobre el tema. Scikit-image también incluye una gran variedad de algoritmos para procesado de imágenes. A partir de ahí, es frecuente de nuevo utilizar sklearn o una biblioteca de deep learning para aprender de las imágenes.
  • SymPy es una biblioteca de cálculo simbólico muy útil en ciertos escenarios relacionados con el aprendizaje automático, al igual que Theano (pero más sencilla).
  • XGBoost es una biblioteca muy popular creada por la comunidad para formular y entrenar modelos de Gradient Boosted Trees. Es más escalable que la que viene en sklearn gracias a su muy ingeniosa implementación. Además, incluye una novedosa regularización que evita muy bien el sobreajuste en la práctica. Se han ganado muchos concursos de Kaggle utilizándola.
  • Scipy incluye un módulo de optimización llamado scipy.optimize que puede ser muy útil cuando formulamos un modelo desde cero, y queremos tener un buen optimizador de la función.
  • StatsModels es una biblioteca estadística y econométrica que utiliza Numpy y Scipy. Personalmente no me gusta mucho, puesto que prefiero utilizar scipy.stats directamente. Pero viene bien explorarla y ver de qué es capaz.
  • Spark (MLlib): Spark tiene una maravillosa API de Python que permite analizar y aplicar aprendizaje automático sobre datasets muy grandes de forma distribuida. La biblioteca MLlib de aprendizaje automático está inspirada en sklearn, hasta el punto que nos encontramos los métodos .fit(), .predict() y .transform() por doquier. La Spark ML Pipeline es también muy similar a las Pipelines de sklearn.
  • Ibis pretende convertirse en un pandas distribuido, exponiendo una API muy similar a la de éste. Utiliza principalmente Impala como motor, y se prevé que en el futuro incluya modelos de machine learning. El ecosistema Blaze pretende hacer algo similar con Dask, pero menos atado al ecosistema Hadoop (de hecho, ya se pueden utilizar algunos modelos de sklearn en distribuido utilizando Dask).

Pero como digo siempre: ¡hay muchas más!

Datahack Logo

Profesor: Julio Antonio Soto